Elephant, neo, quantities, and Spike to Field Coherence

Authors

Dr. Atle Rimehaug | Dr. Nicholas Del Grosso

Download Materials

Setup

Import Libraries

import numpy as np
import quantities as pq
from tqdm import tqdm
import neo
from elephant.sta import spike_field_coherence
from pathlib import Path
import xarray as xr

Utility Functions

Run the cell below, but feel free to ignore this code — it’s some utility functions we’ll use later for convenience.

def calc_wavelet_spectrum(sig, sampling_frequency, wavelet = "cmor1.5-1.0", widths = np.logspace(start=3, stop=10, num=50, base=2)):
    '''
    Calculates the time-resolved power spectrum of a signal using morlet wavelets

    Args:
        sig: signal that the wavelet spectrum is calculated for
        sampling_frequency: sampling frequency of signal
        wavelet: the type of wavelet (bandwidth and ...) to apply in transform
        widths: scales (inverse of frequencies) that the spectrum are calculated for
    Returns:
        freqs: frequencies for which the the spectral values are calculated
        spectrum: complex time-resolved spectrum containing both amplitude and phase over time

    
    '''
    
    sampling_period = 1/sampling_frequency

    spectrum, freqs = pywt.cwt(sig, widths, wavelet, sampling_period=sampling_period)

    return freqs, spectrum

def convert_spikes_to_timestamps(spike_counts):
    '''
    Takes array containing spike counts of units across trials and time and converts the data
    to timestamps of time points of spikes for each unit and trial

    Args:
        spike_counts: array with spike count data
    returns:
        timestamps_all_units: dictionary with timestamps organized by units, trials, time
    
    '''
    timestamps_all_units = dict()
    for iunit in range(spike_counts.shape[0]):
        timestamps_trials = dict()
        for itrial in range(spike_counts.shape[1]):
            timestamps= np.where(spike_counts[iunit,itrial].values)[0]

            timestamps_trials[itrial] = timestamps

        timestamps_all_units[iunit] = timestamps_trials

    return timestamps_all_units

def create_neo_lfp(lfp, sampling_frequency):
    """Put lfp data in time window in neo format"""
    return neo.AnalogSignal(
        lfp.values.T,    # The lfp must be transposed to have time on the first dimension for neo
        units=pq.V, 
        sampling_rate=sampling_frequency*pq.Hz
    )

def get_coherence_all_trials(lfp, timestamps, sampling_frequency_lfp):
    """Loops through all neurons in an area and all trials and calculates the spike field coherence """
    sfc_all_cells = []
    for icell in tqdm(timestamps.keys()):
        sfc_all_trials = []
        for itrial in range(lfp.shape[1]):

            lfp_neo = neo.AnalogSignal(lfp[:,itrial].values.T, units=pq.V, sampling_rate=sampling_frequency_lfp*pq.Hz)

            spiketrain_neo = neo.core.SpikeTrain(timestamps[icell][itrial], units=pq.ms, t_stop = 1.5*pq.s)

            # calculate spike field coherence
            sf_coherence, freqs = spike_field_coherence(lfp_neo, spiketrain_neo)

            sfc_all_trials.append(sf_coherence.T)

        sfc_all_cells.append(np.array(sfc_all_trials))

    return np.array(sfc_all_cells)

class utils:
    calc_wavelet_spectrum = calc_wavelet_spectrum
    convert_spikes_to_timestamps = convert_spikes_to_timestamps
    create_neo_lfp = create_neo_lfp
    get_coherence_all_trials = get_coherence_all_trials

Download Dataset

import owncloud

owncloud.Client.from_public_link('https://uni-bonn.sciebo.de/s/M56DCkpyCJ9FRoi', folder_password="ibots"
).get_file('/', 'dataset_session_754312389.nc')

True

Introduction

In this notebook, you will learn how to use the libraries elephant, neo, and quantities when working with electrophysiological data as well as how to calculate the spike-field coherence between recorded spikes and LFP.

import quantities as pq
import neo

import numpy as np
import matplotlib.pyplot as plt

Section 1: Quantities. Adding units to variables

Checking the units of your final answer in a calculation is a good way to make sure that your calculation is correct. If your answer has the wrong units, then there’s definitely something wrong with the calculation. In Python, we typically don’t have units attached to variables, so we have to keep track of them through a calculation ourselves. However, there are libraries that allow you to attach units to variables, so that it can tracked automatically. One of those libraries is quantities, which we’ll be working with in this section.

Code	Description
`pq.A`	Creates the quantities unit Ampere
`pq.V`	Creates the quantities unit Volt
`pq.Hz`	Creates the quantities unit Hertz
`pq.s`	Creates the quantities unit seconds
`pq.ms`	Creates the quantities unit milliseconds
`pq.`	If you only write “`pq.`”, you will see a list of all the functions and attributes belonging to `pq`. The attributes will be the different units you can create. Scrolling through the list is one way to find the unit you need.
`x = 10*pq.Hz`	Creates the variable `x` containing the value `10` with the unit `Hz` (Hertz) attached.
`x = np.array([1,2,3] * pq.ms)`	Create a variable `x` containing three values with units milliseconds.
`x.units`	Get the units of the variable `x`.
`x.magnitude`	Get the magnitude (value) of the variable `x` without the unit.
`x * y` or `x / y`	Multiplication or division between two variables with units will also multiply or divide the units of those variables.

Exercises

Example: Create a variable named current that has the value 5 and unit Ampere (A). Display the current variable.

The output should be array(5.) * A, which says that it’s an array containing the value 5 multiplied with the unit A (for Ampere).

current = 5*pq.A

current

array(5.) * A

Exercise: Create a variable named potential which contains the value 10 and has the unit Volt. Display the potential variable.

Solution

potential = 10 * pq.V

potential

array(10.) * V

Exercise: Create and display a variable named time_points which contains the three values [10,20,30] and has the unit seconds.

Solution

time_points = [10,20,30] * pq.s

time_points

array([10., 20., 30.]) * s

Exercise: Get the units of the potential variable. The output should be array(1.) * V.

Solution

potential.units

array(1.) * V

Exercise: Get the units of the time variable. The output should be array(1.) * s.

Solution

time_points.units

array(1.) * s

Exercise: Get the values (the magnitude) of the time_points variable (that is, the values of the array without the unit). The output should be array([10., 20., 30.]).

Solution

time_points.magnitude

array([10., 20., 30.])

Exercise: Multiply distance a with distance b. What’s the resulting unit?

a = 10*pq.m
b = 2*pq.m

Solution

a*b

array(20.) * m**2

Exercise: Divide distance a with distance b. What unit does the result have?

a = 10*pq.m
b = 2*pq.m

Solution

a/b

array(5.) * dimensionless

Answer: The result is dimensionless, meaning that it has no unit.

Section 2: Neo: Representing Electrophysiology Data in Python

With Neo, you can store data, metadata, as well as the relationships between different parts of the data together in one Neo object.

Representing time series data

AnalogSignal: Continuous data sampled in regular intervals
IrregularlySampledSignal: Continuous data sampled in irregular intervals

Code	Description
`anasig = neo.AnalogSignal(signal, units = pq.V, sampling_rate, t_start)`	Creates an analog signal neo object containing continuous signal data sampled at regular intervals. Required metadata: physical unit of samples, timestamps of the samples (first timestamp and sampling interval)
`anasig.annotate(quality='excellent', date='2030-05-27)`	Add annotations (additional, optional metadata) to the neo AnalogSignal object.
`spiketrain = neo.SpikeTrain(timestamps, units='ms',t_stop=300*pq.ms)`	Make a neo SpikeTrain object containing the timestamps at which a cell spiked. `t_stop` specifies the end of the time window (the highest possible timestamp value) in units of milliseconds (`pq.ms`).
`neo_object.magnitude`	Get the values (the magnitude of the signal, for example) of the neo object without the unit.
`neo_object.units`	Get the units of neo object.
`np.random.random(size=(nrows,ncols))`	Creates an array containing random values between 0 and 1 with `nrows` number of rows and `ncols` number of columns.

Exercises

Example: Create an AnalogSignal object out of the recorded_potential variable representing a fictional LFP trace from a single channel. The signal object should have units V (Volt), a sampling rate of 1000 Hz, and start at 0 ms.

recorded_potential = [0,1,2,3,4,5,5,4,3,2,1,0]

lfp_trace = neo.AnalogSignal(recorded_potential,
                            units='V', 
                            sampling_rate=1000.*pq.Hz,
                            t_start=0*pq.ms)

lfp_trace

AnalogSignal with 1 channels of length 12; units V; datatype int64
sampling rate: 1000.0 Hz
time: 0.0 ms to 12.0 ms

Exercise: Create an AnalogSignal object of the recorded_potential variable containing fictional LFP traces from 2 recording channels with 50 samples over time. The signal should have units uV (microvolt), a sampling rate of 10000 Hz, and start at 0 ms.

Solution

recorded_potential = np.random.random((50,2))

lfp_trace = neo.AnalogSignal(recorded_potential, 
                          units='uV',
                          sampling_rate=10000*pq.Hz,
                          t_start=0*pq.ms)
lfp_trace

AnalogSignal with 2 channels of length 50; units uV; datatype float64
sampling rate: 10000.0 Hz
time: 0.0 ms to 5.0 ms

Exercise: Annotate your AnalogSignal from either the previous exercise or the example with the recording date (for example, today’s date). Display the variable containing the AnalogSignal to check that the annotations are included.

Solution

# solution
lfp_trace.annotate(date='2024-03-11')
lfp_trace

AnalogSignal with 2 channels of length 50; units uV; datatype float64
annotations: {'date': '2024-03-11'}
sampling rate: 10000.0 Hz
time: 0.0 ms to 5.0 ms

You can also use neo to store spike train data together with metadata in the same object.

Exercise: Create a SpikeTrain object containing the timestamps (the time points at which the neuron fired) [50, 100, 200] of neuron. Give it units ms and have the period start at 0 and stop 300 ms. Display the spike train.

Solution

timestamps = [50, 100, 200]

spike_train = neo.SpikeTrain(timestamps,
                    units='ms',
                    t_start=0,
                    t_stop=300)

spike_train

SpikeTrain containing 3 spikes; units ms; datatype float64 
time: 0.0 ms to 300.0 ms

You can always show the timestamps by typing .magnitude behind the spike train object, like in the cell below (simply run it, you don’t need to add to the code).

spike_train.magnitude

array([ 50., 100., 200.])

Exercise: Create a SpikeTrain object out of the variable timestamps containing the time points at which a neuron fired. Give it units ms and have the period start at 100 and stop 300 ms. Display the spike train. How many spikes does it say that there are?

timestamps = [120,140,160,180,200,220,240]

Solution

spike_train = neo.SpikeTrain(timestamps,
                    units='ms',
                    t_start=100*pq.ms,
                    t_stop=250*pq.ms)

spike_train

SpikeTrain containing 7 spikes; units ms; datatype float64 
time: 100.0 ms to 250.0 ms

Exercise: Get the magnitude of the neo SpikeTrain object.

Solution

spike_train.magnitude

array([120., 140., 160., 180., 200., 220., 240.])

Section 3: Spike-field coherence between LGN spikes and V1 LFP

Coherence is a measure of the linear association between two signals. The spike-field coherence is the association between the spike train of a single neuron and the LFP trace in a single channel and can be used to quantify the relationship between spikes and LFP.

Code	Description
`lfp_neo = utils.create_neo_lfp(lfp_V1[:,itrial], sampling_frequency_lfp)`	Utility function to create an analog signal neo object containing the LFP from a single trial `itrial`.
`spiketrain = neo.core.SpikeTrain(timestamps, t_stop = 1.5*pq.s)`	Make a neo SpikeTrain object containing the timestamps at which a cell spiked. `t_stop` specifies the end of the time window (the highest possible timestamp value) in units of seconds (`pq.s`).
`sf_coherence, freqs = spike_field_coherence(lfp_neo, spiketrain=st)`	Calculate the spike-field coherence (sfc) between a spike train and an LFP trace.
`plt.pcolormesh(x, y, C, cmap = 'plasma', shading = 'gouraud')`	Make a 2D colormap of values in a 2D array (C) against x and y values. The optional parameters `cmap` and `shading` defines, respectively, the colormap and the smoothing applied to the plot.
`plt.plot(x, y)`	Plots the data in `y` against the corresponding values in `x`.
`plt.xlim([some_lower_xlim, some_upper_xlim])`	Limit the x-axis in a plot to be between the numbers `some_lower_xlim` and `some_upper_xlim`.
`transposed_data = data.T`	Transposes the dimensions of the data. An array of dimensions N x M will be transposed to M x N.
`np.nanmean(data, axis = (dim_num))`	Calculate the average of the data across the `dim_num` dimension of the array while ignoring all NaN-values.

Load dataset

Run cells below to load dataset and extract LFP and spike data.

Exercises

loadpath = 'dataset_session_754312389.nc'

dataset_path = Path(loadpath)

dataset = xr.load_dataset(dataset_path)
dataset

<xarray.Dataset> Size: 225MB
Dimensions:                 (channel_depth: 23, trial_nr: 75, time_lfp: 1875,
                             time_csd: 1875, time_whole_rec_lfp: 373124,
                             time_whole_rec_ecp: 373124, unit_id_LGN: 27,
                             time_spikes: 1500, unit_id_V1: 91, unit_id_LM: 13,
                             stimulus_start_times: 75, stimulus_stop_times: 75,
                             channel_id: 23)
Coordinates: (12/13)
  * channel_depth           (channel_depth) int64 184B 0 -40 -80 ... -840 -880
  * trial_nr                (trial_nr) int32 300B 0 1 2 3 4 5 ... 70 71 72 73 74
  * time_lfp                (time_lfp) float64 15kB -1e+03 -999.2 ... 499.2
  * time_csd                (time_csd) float64 15kB -1e+03 -999.2 ... 499.2
  * time_whole_rec_lfp      (time_whole_rec_lfp) float64 3MB 1.286e+06 ... 1....
  * time_whole_rec_ecp      (time_whole_rec_ecp) float64 3MB 1.286e+06 ... 1....
    ...                      ...
  * time_spikes             (time_spikes) float64 12kB -999.5 -998.5 ... 499.5
  * unit_id_V1              (unit_id_V1) int32 364B 951795075 ... 951798053
  * unit_id_LM              (unit_id_LM) int32 52B 951791074 ... 951792163
  * stimulus_start_times    (stimulus_start_times) float64 600B 1.29e+06 ... ...
  * stimulus_stop_times     (stimulus_stop_times) float64 600B 1.29e+06 ... 1...
  * channel_id              (channel_id) int32 92B 850144538 ... 850144362
Data variables:
    lfp_V1                  (channel_depth, trial_nr, time_lfp) float64 26MB ...
    csd_V1                  (channel_depth, trial_nr, time_csd) float64 26MB ...
    lfp_whole_recording_V1  (channel_depth, time_whole_rec_lfp) float64 69MB ...
    ecp_whole_recording_V1  (channel_depth, time_whole_rec_ecp) float64 69MB ...
    spike_counts_LGN        (unit_id_LGN, trial_nr, time_spikes) int16 6MB 0 ...
    spike_counts_V1         (unit_id_V1, trial_nr, time_spikes) int16 20MB 0 ...
    spike_counts_LM         (unit_id_LM, trial_nr, time_spikes) int16 3MB 0 ....
    pupil_width             (trial_nr) float64 600B 39.12 40.53 ... 46.57 44.31
    run_speed               (trial_nr) float64 600B 1.155 1.597 ... 1.251 1.711
Attributes:
    time_unit:                  millisecond
    lfp_unit:                   Volt
    channel_depth_unit:         micrometer
    note_channel_depth:         Measured in distance from electrode closest t...
    sampling_frequency_lfp:     1250
    sampling_frequency_spikes:  1000
    sampling_frequency_unit:    Hz
    stimulus_onset:             1000
    stimulus_offset:            1250

# extract LFP data
lfp_V1 = dataset['lfp_V1']
channels_depth = dataset['channel_depth']
time_trial_lfp = lfp_V1.time_lfp

sampling_frequency_lfp = dataset.sampling_frequency_lfp

# extract spike data
spike_counts_LGN = dataset['spike_counts_LGN']
spike_counts_V1 = dataset['spike_counts_V1']
spike_counts_LM = dataset['spike_counts_LM']

# convert spike counts to timestamps
timestamps_LGN = convert_spikes_to_timestamps(spike_counts_LGN)
timestamps_V1 = convert_spikes_to_timestamps(spike_counts_V1)
timestamps_LM = convert_spikes_to_timestamps(spike_counts_LM)

Example: Compute the spike field coherence between the LFP in trial 0 and the LGN spike train in trial 0 for cell 0. Plot the coherence using plt.pcolormesh.

itrial = 0
lfp_neo = utils.create_neo_lfp(lfp_V1[:,itrial], sampling_frequency_lfp)

icell = 0
timestamps = timestamps_LGN[icell][itrial]
spiketrain_neo = neo.core.SpikeTrain(timestamps, units=pq.ms, t_stop = 1500*pq.ms)

# calculate spike field coherence
sf_coherence, freqs = spike_field_coherence(lfp_neo, spiketrain_neo)

plt.pcolormesh(freqs, channels_depth, sf_coherence.T, cmap = 'inferno', shading = 'gouraud')
plt.xlabel('Frequency (Hz)')
plt.ylabel(r'Depth ($\mu$m)');
plt.title(f'Spike field coherence in trial {itrial} between LGN cell nr. {icell} and V1 LFP')
plt.colorbar()

Exercise: Compute the spike field coherence between the LFP and the LGN spike train in trial at index 5 for cell at index 2. Plot the coherence using plt.pcolormesh. You can use the utility function to create the LFP neo object (found both in the reference table and the example above.)

Note: In this case, the resulting plot will mainly contain noise, so don’t expect to see any structure that you can interpret from the plot.

Solution

itrial = 5
lfp_neo = utils.create_neo_lfp(lfp_V1[:,itrial], sampling_frequency_lfp)

icell = 2
timestamps = timestamps_LGN[icell][itrial]*pq.ms
spiketrain_neo = neo.core.SpikeTrain(timestamps, t_stop = 1500*pq.ms)

# calculate spike field coherence
sf_coherence, freqs = spike_field_coherence(lfp_neo, spiketrain_neo)

plt.pcolormesh(freqs, channels_depth, sf_coherence.T, cmap = 'inferno', shading = 'gouraud')
plt.xlabel('Frequency (Hz)')
plt.ylabel(r'Depth ($\mu$m)');
plt.title(f'Spike field coherence in trial {itrial} between LGN cell nr. {icell} and V1 LFP')
plt.colorbar()

For single trials and single cells it’s often not possible to identify any clear relationships in the spike-field coeherence. If you average over cells, over trials, or both, potential relationships should become more visible.

DEMO: Calculate the spike-field coherence between between LGN spike trains and V1 LFP in trial 0 for all cells. Loop through all cells, calculate the spike-field coherence, and store the result in a list.

Make a colormap plot of the average spike-field coherence across cells with pcolormesh. Are there any frequencies where the coherence seems to be notably higher?

# provided

itrial = 0
lfp_neo = utils.create_neo_lfp(lfp_V1[:,itrial], sampling_frequency_lfp)

sfc_all_cells_LGN = []
for icell in timestamps_LGN.keys():
    spiketrain_neo = neo.core.SpikeTrain(timestamps_LGN[icell][itrial], units=pq.ms, t_stop = 1500*pq.ms)
    # calculate spike field coherence
    sfc, freqs = spike_field_coherence(lfp_neo, spiketrain_neo)

    sfc_all_cells_LGN.append(sfc.T)

sfc_cell_avg = np.nanmean(sfc_all_cells_LGN, axis = 0)

# solution
plt.pcolormesh(freqs, channels_depth, sfc_cell_avg, cmap = 'inferno', shading = 'gouraud')
plt.xlabel('Frequency (Hz)')
plt.ylabel(r'Depth ($\mu$m)');
plt.xlim([0, 200])
plt.colorbar()

Exercise: Compute the average spike field coherence between LGN spikes and V1 LFP (the variable named sfc_all_cells_LGN) across both cells AND channels and plot the result in a 1D against the frequencies. Are the major peak(s) where you would expect if you compare to the 2D colormap plot of coherence across channels above?

Hint: You may have to use np.nanmean to calculate the average, in case there are NaNs in the calculated coherence for some cells.

Solution

sfc_avg = np.nanmean(sfc_all_cells_LGN, axis = (0,1))

plt.plot(freqs, sfc_avg)
plt.xlabel('Frequency (Hz)')
plt.xlim([0,200]);

Example: Calculate the spike-field coherence between between LGN spike trains and V1 LFP for all trials AND for all cells. Compute the average across cells, channels, trials and plot the result against the frequencies.

At which frequencies are the two main peaks in the spike-field coherence between LGN spikes and V1 LFP?

sfc_all_trials_LGN = utils.get_coherence_all_trials(lfp_V1, timestamps_LGN, sampling_frequency_lfp)


sfc_trial_avg_chan_avg_LGN = np.nanmean(sfc_all_trials_LGN, axis = (0, 1, 2))

plt.plot(freqs, sfc_trial_avg_chan_avg_LGN)
plt.xlim([0,200])

plt.xlabel('Frequencies (Hz)')
plt.ylabel('Spectral coherence score (a. u.)')
plt.title('Trial avg. spike-field coherence between LGN spikes\nand V1 LFP. Averaged across all channels');

Let’s check which peaks in the spike-field coherence between LGN spikes and V1 LFP can also be found in the spike-field coherence between spikes in LM or V1 and V1 LFP. If a peak in the spike-field coherence can be observed for all three areas - LGN, V1, and LM - then that’s an indication that it’s an oscillation that permeates the whole brain region. If a peak can only be observed for one area, that indicates a specific relationship between the spikes in that area (at a certain frequency) and the V1 LFP.

Exercise: Calculate the spike-field coherence between LM spike trains and V1 LFP for all trials AND for all cells. Compute the average across cells, channels, trials and plot the result against the frequencies.

Which of the peaks observed in the spike-field coherence between LGN spikes and V1 LFP can also be observed here?

Solution

sfc_all_trials_LM = utils.get_coherence_all_trials(lfp_V1, timestamps_LM, sampling_frequency_lfp)

sfc_trial_avg_chan_avg_LM = np.nanmean(sfc_all_trials_LM, axis = (0, 1, 2))

plt.plot(freqs, sfc_trial_avg_chan_avg_LM)
plt.xlim([0,200])

plt.xlabel('Frequencies (Hz)')
plt.ylabel('Spectral coherence score (a. u.)')
plt.title('Trial avg. spike-field coherence between LM spikes\nand V1 LFP. Averaged across all channels');


  0%|                                                                                                                   | 0/13 [00:00<?, ?it/s]


  8%|████████▏                                                                                                  | 1/13 [00:01<00:18,  1.51s/it]


 15%|████████████████▍                                                                                          | 2/13 [00:02<00:16,  1.49s/it]


 23%|████████████████████████▋                                                                                  | 3/13 [00:04<00:14,  1.46s/it]


 31%|████████████████████████████████▉                                                                          | 4/13 [00:05<00:13,  1.47s/it]


 38%|█████████████████████████████████████████▏                                                                 | 5/13 [00:07<00:11,  1.49s/it]


 46%|█████████████████████████████████████████████████▍                                                         | 6/13 [00:08<00:10,  1.49s/it]


 54%|█████████████████████████████████████████████████████████▌                                                 | 7/13 [00:10<00:08,  1.49s/it]


 62%|█████████████████████████████████████████████████████████████████▊                                         | 8/13 [00:11<00:07,  1.50s/it]


 69%|██████████████████████████████████████████████████████████████████████████                                 | 9/13 [00:13<00:05,  1.49s/it]


 77%|█████████████████████████████████████████████████████████████████████████████████▌                        | 10/13 [00:14<00:04,  1.50s/it]


 85%|█████████████████████████████████████████████████████████████████████████████████████████▋                | 11/13 [00:16<00:03,  1.51s/it]


 92%|█████████████████████████████████████████████████████████████████████████████████████████████████▊        | 12/13 [00:17<00:01,  1.52s/it]


100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████| 13/13 [00:19<00:00,  1.52s/it]


100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████| 13/13 [00:19<00:00,  1.50s/it]

Exercise: Calculate the spike-field coherence between between V1 spike trains and V1 LFP for all trials AND for all cells. Compute the average across cells, channels, trials and plot the result against the frequencies.

Which of the peaks observed in the spike-field coherence between LGN spikes and V1 LFP can also be observed here?

Solution

sfc_all_trials_V1 = utils.get_coherence_all_trials(lfp_V1, timestamps_V1, sampling_frequency_lfp)

sfc_trial_avg_chan_avg_V1 = np.nanmean(sfc_all_trials_V1, axis = (0, 1, 2))

plt.plot(freqs, sfc_trial_avg_chan_avg_V1)
plt.xlim([0,200])

plt.xlabel('Frequencies (Hz)')
plt.ylabel('Spectral coherence score (a. u.)')
plt.title('Trial avg. spike-field coherence between V1 spikes\nand V1 LFP. Averaged across all channels');


  0%|                                                                                                                   | 0/91 [00:00<?, ?it/s]


  1%|█▏                                                                                                         | 1/91 [00:01<02:14,  1.49s/it]


  2%|██▎                                                                                                        | 2/91 [00:03<02:13,  1.50s/it]


  3%|███▌                                                                                                       | 3/91 [00:04<02:12,  1.50s/it]


  4%|████▋                                                                                                      | 4/91 [00:05<02:10,  1.50s/it]


  5%|█████▉                                                                                                     | 5/91 [00:07<02:10,  1.52s/it]


  7%|███████                                                                                                    | 6/91 [00:09<02:10,  1.53s/it]


  8%|████████▏                                                                                                  | 7/91 [00:10<02:09,  1.54s/it]


  9%|█████████▍                                                                                                 | 8/91 [00:12<02:07,  1.53s/it]


 10%|██████████▌                                                                                                | 9/91 [00:13<02:06,  1.54s/it]


 11%|███████████▋                                                                                              | 10/91 [00:15<02:05,  1.55s/it]


 12%|████████████▊                                                                                             | 11/91 [00:16<02:04,  1.56s/it]


 13%|█████████████▉                                                                                            | 12/91 [00:18<02:03,  1.56s/it]


 14%|███████████████▏                                                                                          | 13/91 [00:20<02:02,  1.57s/it]


 15%|████████████████▎                                                                                         | 14/91 [00:21<02:00,  1.57s/it]


 16%|█████████████████▍                                                                                        | 15/91 [00:23<01:59,  1.57s/it]


 18%|██████████████████▋                                                                                       | 16/91 [00:24<01:57,  1.57s/it]


 19%|███████████████████▊                                                                                      | 17/91 [00:26<01:54,  1.55s/it]


 20%|████████████████████▉                                                                                     | 18/91 [00:27<01:53,  1.55s/it]


 21%|██████████████████████▏                                                                                   | 19/91 [00:29<01:50,  1.54s/it]


 22%|███████████████████████▎                                                                                  | 20/91 [00:30<01:49,  1.54s/it]


 23%|████████████████████████▍                                                                                 | 21/91 [00:32<01:49,  1.56s/it]


 24%|█████████████████████████▋                                                                                | 22/91 [00:34<01:47,  1.55s/it]


 25%|██████████████████████████▊                                                                               | 23/91 [00:35<01:45,  1.55s/it]


 26%|███████████████████████████▉                                                                              | 24/91 [00:37<01:42,  1.52s/it]


 27%|█████████████████████████████                                                                             | 25/91 [00:38<01:40,  1.52s/it]


 29%|██████████████████████████████▎                                                                           | 26/91 [00:40<01:39,  1.53s/it]


 30%|███████████████████████████████▍                                                                          | 27/91 [00:41<01:36,  1.51s/it]


 31%|████████████████████████████████▌                                                                         | 28/91 [00:42<01:33,  1.49s/it]


 32%|█████████████████████████████████▊                                                                        | 29/91 [00:44<01:32,  1.49s/it]


 33%|██████████████████████████████████▉                                                                       | 30/91 [00:45<01:30,  1.48s/it]


 34%|████████████████████████████████████                                                                      | 31/91 [00:47<01:27,  1.47s/it]


 35%|█████████████████████████████████████▎                                                                    | 32/91 [00:48<01:25,  1.46s/it]


 36%|██████████████████████████████████████▍                                                                   | 33/91 [00:50<01:25,  1.47s/it]


 37%|███████████████████████████████████████▌                                                                  | 34/91 [00:51<01:23,  1.47s/it]


 38%|████████████████████████████████████████▊                                                                 | 35/91 [00:53<01:22,  1.48s/it]


 40%|█████████████████████████████████████████▉                                                                | 36/91 [00:54<01:21,  1.48s/it]


 41%|███████████████████████████████████████████                                                               | 37/91 [00:56<01:20,  1.49s/it]


 42%|████████████████████████████████████████████▎                                                             | 38/91 [00:57<01:19,  1.49s/it]


 43%|█████████████████████████████████████████████▍                                                            | 39/91 [00:59<01:17,  1.48s/it]


 44%|██████████████████████████████████████████████▌                                                           | 40/91 [01:00<01:15,  1.49s/it]


 45%|███████████████████████████████████████████████▊                                                          | 41/91 [01:02<01:14,  1.48s/it]


 46%|████████████████████████████████████████████████▉                                                         | 42/91 [01:03<01:12,  1.47s/it]


 47%|██████████████████████████████████████████████████                                                        | 43/91 [01:05<01:10,  1.46s/it]


 48%|███████████████████████████████████████████████████▎                                                      | 44/91 [01:06<01:09,  1.48s/it]


 49%|████████████████████████████████████████████████████▍                                                     | 45/91 [01:08<01:09,  1.50s/it]


 51%|█████████████████████████████████████████████████████▌                                                    | 46/91 [01:09<01:07,  1.51s/it]


 52%|██████████████████████████████████████████████████████▋                                                   | 47/91 [01:11<01:07,  1.53s/it]


 53%|███████████████████████████████████████████████████████▉                                                  | 48/91 [01:12<01:06,  1.54s/it]


 54%|█████████████████████████████████████████████████████████                                                 | 49/91 [01:14<01:04,  1.54s/it]


 55%|██████████████████████████████████████████████████████████▏                                               | 50/91 [01:15<01:03,  1.54s/it]


 56%|███████████████████████████████████████████████████████████▍                                              | 51/91 [01:17<01:02,  1.56s/it]


 57%|████████████████████████████████████████████████████████████▌                                             | 52/91 [01:19<01:00,  1.56s/it]


 58%|█████████████████████████████████████████████████████████████▋                                            | 53/91 [01:20<00:59,  1.56s/it]


 59%|██████████████████████████████████████████████████████████████▉                                           | 54/91 [01:22<00:57,  1.55s/it]


 60%|████████████████████████████████████████████████████████████████                                          | 55/91 [01:23<00:55,  1.55s/it]


 62%|█████████████████████████████████████████████████████████████████▏                                        | 56/91 [01:25<00:54,  1.56s/it]


 63%|██████████████████████████████████████████████████████████████████▍                                       | 57/91 [01:26<00:52,  1.55s/it]


 64%|███████████████████████████████████████████████████████████████████▌                                      | 58/91 [01:28<00:51,  1.56s/it]


 65%|████████████████████████████████████████████████████████████████████▋                                     | 59/91 [01:29<00:49,  1.55s/it]


 66%|█████████████████████████████████████████████████████████████████████▉                                    | 60/91 [01:31<00:48,  1.56s/it]


 67%|███████████████████████████████████████████████████████████████████████                                   | 61/91 [01:33<00:46,  1.56s/it]


 68%|████████████████████████████████████████████████████████████████████████▏                                 | 62/91 [01:34<00:45,  1.56s/it]


 69%|█████████████████████████████████████████████████████████████████████████▍                                | 63/91 [01:36<00:43,  1.55s/it]


 70%|██████████████████████████████████████████████████████████████████████████▌                               | 64/91 [01:37<00:41,  1.53s/it]


 71%|███████████████████████████████████████████████████████████████████████████▋                              | 65/91 [01:39<00:39,  1.52s/it]


 73%|████████████████████████████████████████████████████████████████████████████▉                             | 66/91 [01:40<00:37,  1.51s/it]


 74%|██████████████████████████████████████████████████████████████████████████████                            | 67/91 [01:42<00:36,  1.50s/it]


 75%|███████████████████████████████████████████████████████████████████████████████▏                          | 68/91 [01:43<00:34,  1.51s/it]


 76%|████████████████████████████████████████████████████████████████████████████████▎                         | 69/91 [01:45<00:33,  1.50s/it]


 77%|█████████████████████████████████████████████████████████████████████████████████▌                        | 70/91 [01:46<00:31,  1.50s/it]


 78%|██████████████████████████████████████████████████████████████████████████████████▋                       | 71/91 [01:48<00:29,  1.49s/it]


 79%|███████████████████████████████████████████████████████████████████████████████████▊                      | 72/91 [01:49<00:28,  1.50s/it]


 80%|█████████████████████████████████████████████████████████████████████████████████████                     | 73/91 [01:51<00:26,  1.50s/it]


 81%|██████████████████████████████████████████████████████████████████████████████████████▏                   | 74/91 [01:52<00:25,  1.51s/it]


 82%|███████████████████████████████████████████████████████████████████████████████████████▎                  | 75/91 [01:54<00:24,  1.52s/it]


 84%|████████████████████████████████████████████████████████████████████████████████████████▌                 | 76/91 [01:55<00:22,  1.53s/it]


 85%|█████████████████████████████████████████████████████████████████████████████████████████▋                | 77/91 [01:57<00:21,  1.53s/it]


 86%|██████████████████████████████████████████████████████████████████████████████████████████▊               | 78/91 [01:58<00:19,  1.52s/it]


 87%|████████████████████████████████████████████████████████████████████████████████████████████              | 79/91 [02:00<00:18,  1.51s/it]


 88%|█████████████████████████████████████████████████████████████████████████████████████████████▏            | 80/91 [02:01<00:16,  1.50s/it]


 89%|██████████████████████████████████████████████████████████████████████████████████████████████▎           | 81/91 [02:03<00:14,  1.49s/it]


 90%|███████████████████████████████████████████████████████████████████████████████████████████████▌          | 82/91 [02:04<00:13,  1.50s/it]


 91%|████████████████████████████████████████████████████████████████████████████████████████████████▋         | 83/91 [02:06<00:12,  1.51s/it]


 92%|█████████████████████████████████████████████████████████████████████████████████████████████████▊        | 84/91 [02:07<00:10,  1.52s/it]


 93%|███████████████████████████████████████████████████████████████████████████████████████████████████       | 85/91 [02:09<00:09,  1.51s/it]


 95%|████████████████████████████████████████████████████████████████████████████████████████████████████▏     | 86/91 [02:10<00:07,  1.52s/it]


 96%|█████████████████████████████████████████████████████████████████████████████████████████████████████▎    | 87/91 [02:12<00:06,  1.53s/it]


 97%|██████████████████████████████████████████████████████████████████████████████████████████████████████▌   | 88/91 [02:13<00:04,  1.53s/it]


 98%|███████████████████████████████████████████████████████████████████████████████████████████████████████▋  | 89/91 [02:15<00:03,  1.54s/it]


 99%|████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 90/91 [02:16<00:01,  1.52s/it]


100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████| 91/91 [02:18<00:00,  1.51s/it]


100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████| 91/91 [02:18<00:00,  1.52s/it]

Section 4: Elephant for analyzing spike data

Elephant is a library in Python made specifically for analyzing electrophysiological data. In this section, you’ll explore some of its functions to calculate spike statistics.

Code	Description
`spike_train = homogeneous_poisson_process(rate=some_ratepq.Hz, t_start=0.pq.s, t_stop=10.*pq.s)`	Generates a poisson process spiketrain firing at a rate `some_rate` for 10 seconds.
`spike_train2 = homogeneous_gamma_process(a=3, b=10pq.Hz, t_start=0.pq.s, t_stop=10.*pq.s)`	Generates a gamma process spiketrain firing at a rate `some_rate` for 10 seconds.
`stats.mean_firing_rate(spike_train)`	Calculates the average firing rate of the spiketrain.
`isi = stats.isi(spike_train)`	Calculates the inter-spike interval (ISI) of the spiketrain.
`stats.cv2(isi)`	Calculates the spike interval variability (CV2 score) for a spike train.
`plt.eventplot([spike_train1.magnitude, spike_train2.magnitude])`	Creates a raster plot for two spike trains.
`stats.time_histogram(spike_trains=spike_train1, bin_size=1*pq.ms)`	Computes the spike counts per unit time (for example per millisecond) from timestamps.

Exercises

from elephant.spike_train_generation import homogeneous_poisson_process, homogeneous_gamma_process
from elephant import statistics as stats

Generate spike data

Run the cell below to generate three spike trains using Elephant. Two spike trains are generated by a poisson process and the last one is generated by a gamma process.

np.random.seed(2)  
spike_trainA = homogeneous_poisson_process(rate=3*pq.Hz, t_start=0.*pq.s, t_stop=20.*pq.s)
spike_trainB = homogeneous_poisson_process(rate=6*pq.Hz, t_start=0.*pq.s, t_stop=20.*pq.s)
spike_trainC = homogeneous_gamma_process(a=3, b=6*pq.Hz, t_start=0.*pq.s, t_stop=20.*pq.s)

spike_trainA

SpikeTrain containing 68 spikes; units s; datatype float64 
time: 0.0 s to 20.0 s

Run the cell below to create a raster plot of the three spike trains.

plt.figure(figsize=(8, 3))
plt.eventplot([spike_trainA.magnitude, spike_trainB.magnitude, spike_trainC.magnitude])
plt.xlabel('Time (s)')
plt.yticks([0, 1, 2], labels=["Poisson process spike train A", "Poisson process spike train B", "Gamma process spike train C"]);

Example: Calculate average firing rate of spike train A.

firing_rate = stats.mean_firing_rate(spike_trainA)
firing_rate

array(3.4) * 1/s

Exercise: Calculate the average firing rate of spike train B.

Solution

firing_rate = stats.mean_firing_rate(spike_trainB)
firing_rate

array(6.55) * 1/s

You can calculate the inter-spike intervals (ISI) for a spike train where you get the unit included (it should say * s at the end of the array containing the ISIs when you display it).

Exercise: Calculate the ISI for spike train A using elephant’s function stats.isi.

Solution

isi_spike_trainA = stats.isi(spike_trainA)

isi_spike_trainA

array([0.00875608, 0.26591931, 0.19050011, 0.18178717, 0.13365914, ...,
       0.06179894, 0.11558582, 0.24749264, 0.1470088 , 0.01558505]) * s

Exercise: Calculate the ISI for spike train B.

Solution

isi_spike_trainB = stats.isi(spike_trainB)

isi_spike_trainB

array([0.07202087, 0.1334472 , 0.60143601, 0.01994301, 0.0626329 , ...,
       0.28215394, 0.02347772, 0.07093529, 0.03897143, 0.33671316]) * s

Exercise: Calculate the ISI for spike train C.

Solution

isi_spike_trainC = stats.isi(spike_trainC)

isi_spike_trainC

array([0.78341479, 0.12469095, 0.3933329 , 0.68786487, 0.2856468 , ...,
       0.09028663, 0.70179214, 0.27786033, 0.31526487, 0.26520429]) * s

The interspike intervals can be used to determine whether a spiking neuron is in a bursting, tonic, or poisson firing mode. A bursting neuron fires very often in a short time window but is otherwise mostly silent. A tonic firing mode means that the neuron is firing at highly regular intervals. A poisson firing mode means that the neuron is firing at random and is neither bursting nor in tonic firing mode.

A metric to calculate the the spike interval variability called the CV2 score (stats.CV2) is computed from the interspike interval of a neuron. If the CV2 score is 1 (or close to 1), the neuron is firing in a poisson firing mode. If it’s > 1, then the neuron is bursting. If it’s 0 (or close to 0), then the neuron is in tonic firing mode.

Example: Calculate the CV2 of the inter-spike interval (ISI) spike train A.

cv2_spike_trainA = stats.cv2(isi_spike_trainA)

cv2_spike_trainA

0.9653690108731857

Exercise: Calculate the CV2 of the inter-spike interval (ISI) spike train B.

Solution

cv2_spike_trainB = stats.cv2(isi_spike_trainB)

cv2_spike_trainB

0.9636532783667118

Exercise: Calculate the CV2 of the inter-spike interval (ISI) spike train C.

Which spike trains have a CV2 score closest to 1? Is it the ones you would have expected considering which spike trains were generated with a poisson process?

Solution

cv2_spike_trainC = stats.cv2(isi_spike_trainC)

cv2_spike_trainC

0.7245379330293805