Courses

Combining Firing Rate and CSD Analysis

Authors

Dr. Atle Rimehaug | Dr. Nicholas Del Grosso

Setup

Import Libraries

import numpy as np
import neo
import quantities as pq
from elephant.current_source_density import estimate_csd
import xarray as xr
from pathlib import Path

Utility Functions

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

def add_um_units(channels_depth_vals):
    """Add units to coordinates of channels on probe"""
    return np.array([np.abs(channels_depth_vals)]).T*pq.um

def create_neo_lfp(lfp_values, sampling_frequency):
    """Put lfp data in time window in neo format"""
    return neo.AnalogSignal(lfp_values, units=pq.V, sampling_rate=sampling_frequency*pq.Hz)

'''def get_data_time_window(time_window_start, time_window_end, data, times_sampled):
    """get lfp and spike data in selected time window (unsure if going to use this because currently need to take different time dimension of 
       csd/lfp and spike data into account)"""
    # get csd and time points of lfp samples in time window
    mask_window = (times_sampled > time_window_start) & (times_sampled < time_window_end)
    time_window = times_sampled[mask_window]
    data_window = data[:, mask_window]

    # get population firing rate and time points of spike data samples in time window
    mask_window_spikes = (time_spikes > time_window_start) & (time_spikes < time_window_end)
    time_window_spikes = time_spikes[mask_window_spikes]
    pop_firing_rate_LGN_window = pop_firing_rate_LGN_trial_avg[mask_window_spikes]'''

class utils:
    add_um_units = add_um_units
    create_neo_lfp = create_neo_lfp

Download Dataset

import owncloud

owncloud.Client.from_public_link('https://uni-bonn.sciebo.de/s/rY4YCPXmkpespKG').get_file('/', 'dataset_session_754312389.nc')

True

Load dataset

Run cells below to load dataset and extract relevant data

Introduction

LFP reflects currents caused by synaptic input. Spikes in a presynaptic population activate the synaptic currents in the postsynaptic population. Therefore, the LFP, and by extension the CSD, reflects the spikes in a presynaptic population. In this notebook, you will learn how to calculate and plot firing rates and combine firing rate analysis with CSD analysis to observe this relationship between spikes and LFP/CSD.

Section 1: Calculate and plot firing rates from spike data

Code	Description
`np.arange(start, stop, step)`	Makes an array with all the integers between two values, with a given spacing
`np.mean(data, axis = (dim_num))` or `data.mean(axis = dim_num)`	Calculate the average of the data across the `dim_num` dimension of the array.
`np.mean(data, axis = (dim_num1, dim_num2))` or `data.mean(axis = (dim_num1, dim_num2)`	Calculate the average of the data across both the `dim_num1` and the `dim_num2` dimensions of the array.
`data.mean('trial_nr')`	If `data` is an `xarray DataArray`, the coordinate names (e.g. `'trial_nr'`) can be used to specify which dimension the average should be calculated for.
`plt.plot(x, y)`	Plot y values against x values.
`second_row = data[1, :]`	Get data on second row in 2D data array.
`third_column = data[:, 2]`	Get data on third column in 2D data array.
`element_second_row_third_column = data[1, 2]`	Get the single element on second row and third column in 2D data array.

Load Dataset

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

Exercises

dataset_path = Path('dataset_session_754312389.nc')

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 data
spike_counts_LGN = dataset['spike_counts_LGN']
spike_counts_V1 = dataset['spike_counts_V1']
spike_counts_LM = dataset['spike_counts_LM']

lfp_V1 = dataset['lfp_V1']
lfp_trial_avg = np.mean(lfp_V1, axis=1)
channels_depth = dataset['channel_depth']
sampling_frequency_lfp = dataset.sampling_frequency_lfp
lfp_V1

<xarray.DataArray 'lfp_V1' (channel_depth: 23, trial_nr: 75, time_lfp: 1875)> Size: 26MB array([[[-8.51569421e-06, -1.15631923e-05, -3.76924291e-07, …, -2.17654367e-06, 1.68119712e-07, -6.94200500e-06], [ 2.81771330e-08, 6.29403676e-06, 9.98029865e-06, …, 1.01304291e-05, 5.93295341e-06, 1.72969475e-05], [-1.27271427e-06, -3.68783112e-06, -4.46340366e-06, …, 3.65666426e-06, 2.96670978e-06, 1.24568597e-06], …, [ 4.70565253e-06, 1.94317598e-06, 1.20455675e-06, …, 4.07503496e-06, 7.65054254e-06, 3.91673831e-06], [ 8.38144661e-07, 2.28123812e-06, 3.70896321e-06, …, 1.34013793e-05, 1.15462596e-05, 8.42411093e-06], [ 7.32479652e-06, 1.07711181e-06, 7.00885884e-06, …, 1.92065863e-05, 1.92065863e-05, 1.92065863e-05]],

   [[ 3.91952787e-06,  1.70056905e-06,  5.97275690e-06, ...,
      4.06974274e-06,  4.04094749e-06, -4.68055730e-06],
    [-4.10311367e-06,  5.00856269e-06,  6.20500275e-06, ...,
      8.75868311e-06,  5.69604080e-06,  8.57744575e-06],
    [ 1.09646506e-06,  1.62070262e-06,  1.10953040e-07, ...,
      1.01990671e-05, -8.36754580e-06, -1.43387447e-05],

… [-6.46738552e-06, -1.19701456e-05, -2.41882630e-05, …, -2.31175915e-05, -2.37374262e-05, -7.53926446e-06], [ 2.91850049e-05, 2.81713398e-05, 5.17618757e-05, …, 4.44764295e-05, 4.08778688e-05, 3.83233192e-05], [-1.80642354e-05, -1.06068608e-05, -3.71177368e-06, …, 1.48063260e-05, 1.48063260e-05, 1.48063260e-05]],

   [[ 4.42062699e-05,  5.99208645e-05,  7.60618853e-05, ...,
      9.85508013e-05,  8.87923195e-05,  8.02998878e-05],
    [-6.30930699e-06,  2.14882295e-05,  5.40961203e-06, ...,
     -6.18517313e-05, -8.49961217e-05, -7.07025406e-05],
    [ 1.32347041e-04,  9.66830285e-05,  9.77152893e-05, ...,
     -3.68393584e-05, -2.08663734e-05, -1.27901811e-05],
    ...,
    [ 7.59530363e-07, -9.00525049e-06, -1.95382678e-05, ...,
     -2.78150664e-05, -1.46349189e-05,  3.65395950e-06],
    [ 3.22381374e-05,  3.90163375e-05,  6.39709157e-05, ...,
      5.17422904e-05,  4.94353600e-05,  3.90770304e-05],
    [-2.50190016e-05, -1.82145705e-05, -1.97245214e-05, ...,
      6.13148062e-06,  6.13148062e-06,  6.13148062e-06]]])

Coordinates:

channel_depth (channel_depth) int64 184B 0 -40 -80 -120 … -800 -840 -880
trial_nr (trial_nr) int32 300B 0 1 2 3 4 5 6 … 68 69 70 71 72 73 74
time_lfp (time_lfp) float64 15kB -1e+03 -999.2 -998.4 … 498.4 499.2

from scipy.ndimage import gaussian_filter1d

import matplotlib.pyplot as plt

Example: Plot the spike count (also often called spike train) of the LGN cell at index 0 in trial index 2.

icell = 0
itrial = 2

time_spikes = spike_counts_LGN.time_spikes
plt.plot(time_spikes, spike_counts_LGN[icell, itrial])
plt.xlabel('Time from stimulus onset (ms)');

Exercise: Plot the spike count of a V1 cell at index 0 in trial index 2.

Solution

icell = 0
itrial = 2

time_spikes = spike_counts_V1.time_spikes
plt.plot(time_spikes, spike_counts_V1[icell, itrial])
plt.xlabel('Time from stimulus onset (ms)');

Exercise: Plot the spike count of an LM cell at index 6 in trial index 50.

Solution

icell = 6
itrial = 50

time_spikes = spike_counts_LM.time_spikes
plt.plot(time_spikes, spike_counts_LM[icell, itrial])
plt.xlabel('Time from stimulus onset (ms)');

Exercise: Calculate the trial averaged firing rate of all LGN cells and plot the firing rate of cell nr. 0.

Solution

icell = 0

firing_rates_LGN = spike_counts_LGN.mean('trial_nr')

plt.plot(time_spikes, firing_rates_LGN[icell]);

Exercise: Calculate and plot the trial averaged firing rates of cell at index 10 in V1.

Does this cell appear to be responding most strongly to the onset and offset of the stimulus? (Onset is at 0 ms and offset is at 250 ms.)

Solution

icell = 10

firing_rates_V1 = spike_counts_V1.mean('trial_nr')

plt.plot(time_spikes, firing_rates_V1[icell]);

Exercise: Calculate and plot the trial averaged firing rates for the cell at index 6 in LM.

Does this cell respond to stimulus onset or offset or both?

Solution

icell = 6

firing_rates_LM = spike_counts_LM.mean('trial_nr')

plt.plot(time_spikes, firing_rates_LM[icell]);

Exercise: Calculate and plot the population firing rates in trial 0 in LGN.

Hint: Useful code to specify which dimension to compute the average for: spike_counts_LGN.mean('unit_id_LGN')

Solution

pop_firing_rate_LGN = spike_counts_LGN.mean('unit_id_LGN')

itrial = 0
plt.plot(time_spikes, pop_firing_rate_LGN[itrial]);

Exercise: Calculate and plot trial averaged population firing rate in LGN. In other words, calculate the average across both all trials and all cells in LGN.

Hint: To compute the average across multiple dimensions of an xarray DataArray, include a list of the coordinate/dimension names to average over, e.g. data.mean(['coord_nameA', 'coord_nameB']).

Is the timing of the firing rate peaks occuring where you would expect?

Solution

pop_firing_rate_LGN = spike_counts_LGN.mean(['trial_nr', 'unit_id_LGN'])

plt.plot(time_spikes, pop_firing_rate_LGN);

Exercise: Calculate and plot trial averaged population firing rate in V1.

Solution

pop_firing_rate_V1 = spike_counts_V1.mean(['trial_nr', 'unit_id_V1'])

plt.plot(time_spikes, pop_firing_rate_V1)

Section 2: Make Figure with Subplots where Firing Rates and CSD Plots are Shown Together.

Code	Description
`plt.plot(x, y)`	Plot y values against x values.
`plt.pcolormesh(x, y, C, cmap='RdBu_r')`	Make a 2D colormap of values in a 2D array (C) against x and y values with the colormap set with `cmap`.
`fig, ax = plt.subplots(nrows, ncols)`	Tells python to create a grid of subplots with `nrows` number of rows and `ncols` number of columns. `ax` can be indexed to refer to a specific subplot, f.ex. the subplot on the first row and second column would be `ax[0,1]`.
`gaussian_filter1d(data, sigma = filter_width)`	Apply a 1D gaussian filter that smoothes the data to reduce the prominence of noise in a plot. `sigma = filter_width` designates the width of the filter.

Exercises

from elephant.current_source_density import estimate_csd

Exercise: Run the cell below to calculate the CSD and the population firing rates. (No need to add or change anything in the code.)

# provided

lfp_neo = utils.create_neo_lfp(lfp_values=lfp_trial_avg.values.T, sampling_frequency=sampling_frequency_lfp) # The lfp must be transposed to have time on the first dimension for neo
coords = utils.add_um_units(channels_depth_vals=channels_depth.values)
csd_neo = estimate_csd(lfp_neo, coordinates=coords, method='DeltaiCSD', diam=400*pq.um).T # The csd must be transposed to get channels on first dimension for the plot

pop_firing_rate_LGN_trial_avg = spike_counts_LGN.mean(['trial_nr', 'unit_id_LGN'])*1E3 # multiply with 1E3 to scale to Hz

# get sampling time points of LFP and spike data
time_lfp = lfp_trial_avg.time_lfp
time_spikes = spike_counts_LGN.time_spikes
csd_neo

discrete filter coefficients: 
b = [ 0.607 1.000 0.607 ],                
a = [ 2.213 ]

AnalogSignal with 1875 channels of length 23; units A/m**2; datatype float64
annotations: {'x_coords': array([[  0.],
         [ 40.],
         [ 80.],
         [120.],
         [160.],
         [200.],
         [240.],
         [280.],
         [320.],
         [360.],
         [400.],
         [440.],
         [480.],
         [520.],
         [560.],
         [600.],
         [640.],
         [680.],
         [720.],
         [760.],
         [800.],
         [840.],
         [880.]]) * um}
sampling rate: 1250.0 Hz
time: 0.0 s to 0.0184 s

Example: Make a figure with one subplot containing the V1 CSD and one subplot containing the trial averaged LGN firing rate from 500 ms before stimulus onset to 500 ms after stimulus onset (both the firing rate and the CSD in this time window are extracted for you if you run the cell below).

# select time window

time_window_start = -500
time_window_end = 500

# get csd and time points of lfp samples in time window
mask_window_lfp = (time_lfp > time_window_start) & (time_lfp < time_window_end)
time_window_lfp = time_lfp[mask_window_lfp]
csd_window = csd_neo[:, mask_window_lfp]

# get population firing rate and time points of spike data samples in time window
mask_window_spikes = (time_spikes > time_window_start) & (time_spikes < time_window_end)
time_window_spikes = time_spikes[mask_window_spikes]
pop_firing_rate_LGN_window = pop_firing_rate_LGN_trial_avg[mask_window_spikes]

# Plot CSD and LGN population firing rates together with subplots

fig, ax = plt.subplots(nrows = 2)

ax[0].pcolormesh(time_window_lfp, channels_depth, csd_window, cmap = 'RdBu_r')
ax[0].set_title('CSD in V1')

ax[1].plot(time_window_spikes, pop_firing_rate_LGN_window)
ax[1].set_title('Firing rate LGN')

fig.tight_layout()

Exercise: Make a figure with one subplot containing the V1 CSD and one subplot containing the trial averaged LGN firing rate for the time window from 0 ms before stimulus onset to 100 ms after (both are calculated for you in the cell below).

Do you notice anything when you compare the timing of the first peak in the LGN firing rate and the onset of the first sink in the V1 CSD?

Solution

# select time window

'''Set time_window_start and time_window_end here'''
time_window_start = 0
time_window_end = 100

# get csd and time points of lfp samples in time window
mask_window_lfp = (time_lfp > time_window_start) & (time_lfp < time_window_end)
time_window_lfp = time_lfp[mask_window_lfp]
csd_window = csd_neo[:, mask_window_lfp]

# get population firing rate and time points of spike data samples in time window
mask_window_spikes = (time_spikes > time_window_start) & (time_spikes < time_window_end)
time_window_spikes = time_spikes[mask_window_spikes]
pop_firing_rate_LGN_window = pop_firing_rate_LGN_trial_avg[mask_window_spikes]

# solution

# Plot CSD and LGN population firing rates together with subplots

fig, ax = plt.subplots(nrows = 2, figsize = (3,6))

ax[0].pcolormesh(time_window_lfp, channels_depth, csd_window, cmap = 'RdBu_r')
ax[0].set_title('CSD in V1')

ax[1].plot(time_window_spikes, pop_firing_rate_LGN_window)
ax[1].set_title('Firing rate LGN')

fig.tight_layout()

Due to noise, it’s not always that easy to determine what the initial peak is when looking at such firing rate plots. Smoothing with filters can help reduce the noise in the plot and make the onset of peaks more visible.

Exercise: Apply a gaussian filter with sigma = 1 to the firing rate for LGN. Then make the same figure with subplots for the trial averaged population firing rate in LGN and the CSD for the first 100 ms after stimulus onset.

Solution

# Plot CSD and LGN population firing rates together with subplots

fig, ax = plt.subplots(nrows = 2, figsize = (3,6))

ax[0].pcolormesh(time_window_lfp, channels_depth, csd_window, cmap = 'RdBu_r')
ax[0].set_title('CSD in V1')

ax[1].plot(time_window_spikes, gaussian_filter1d(pop_firing_rate_LGN_window, sigma=1))
ax[1].set_title('Firing rate LGN')

fig.tight_layout()

Exercise: Make a figure with the CSD of V1 for the first 100 ms after stimulus onset in one subplot and the trial-averaged firing rates of the higher visual area LM for the same time period in another subplot.

Hint: You need to first calculate the average across all cells and all trials to get the population firing rate in LM, then apply the mask for the 100 ms time window to slice out the firing rate in this time period. The code to do both these operations to LGN spikes can be found above.

Without smoothing, do you notice anything when you compare the timing of the first peak in the LM firing rate and the onset of the later sinks in the V1 CSD?

Solution

# compute trial- and population averaged firing rate of LM
pop_firing_rate_LM_trial_avg = spike_counts_LM.mean(['trial_nr', 'unit_id_LM'])*1E3 # multiply with 1E3 to scale to Hz

# get population firing in time window
pop_firing_rate_LM_window = pop_firing_rate_LM_trial_avg[mask_window_spikes]

fig, ax = plt.subplots(nrows = 2, figsize = (3,6))

ax[0].pcolormesh(time_window_lfp, channels_depth, csd_window, cmap = 'RdBu_r')
ax[0].set_title('CSD in V1')

ax[1].plot(time_window_spikes, pop_firing_rate_LM_window)
ax[1].set_title('Firing rate LM')

fig.tight_layout()

Exercise: Apply a gaussian filter with sigma = 1 to the firing rate for LM. Then make the same figure with subplots for the trial averaged population firing rate in LGN and the CSD for the first 100 ms after stimulus onset.

Solution

fig, ax = plt.subplots(nrows = 2, figsize = (3,6))

ax[0].pcolormesh(time_window_lfp, channels_depth, csd_window, cmap = 'RdBu_r')
ax[0].set_title('CSD in V1')

ax[1].plot(time_window_spikes, gaussian_filter1d(pop_firing_rate_LM_window, sigma=1))
ax[1].set_title('Firing rate LM')

fig.tight_layout()

Exercise: Make a figure with the CSD of V1 for the first 100 ms after stimulus onset in one subplot and the trial-averaged firing rates of the V1 for the same time period in another subplot. Apply smoothing with sigma = 1.

Solution

# compute trial- and population averaged firing rate of LM
pop_firing_rate_V1_trial_avg = np.mean(spike_counts_V1, axis=(0,1))*1E3 # multiply with 1E3 to scale to Hz

# get population firing rate of LM in time window
pop_firing_rate_V1_window = pop_firing_rate_V1_trial_avg[mask_window_spikes]

fig, ax = plt.subplots(nrows = 2, figsize = (3,6))

ax[0].pcolormesh(time_window_lfp, channels_depth, csd_window, cmap = 'RdBu_r')
ax[0].set_title('CSD in V1')

ax[1].plot(time_window_spikes, gaussian_filter1d(pop_firing_rate_V1_window, sigma=1))
ax[1].set_title('Firing rate V1')

fig.tight_layout()

Section 3: Add Axes to Plots

An alternative to creating a figure with subplots is to plot the CSD and the firing rates together in the same plot. To do that, we need to add a third axis to the plot.

Code	Description
`plt.plot(x, y)`	Plot y values against x values.
`plt.pcolormesh(x, y, C, cmap='RdBu_r')`	Make a 2D colormap of values in a 2D array (C) against x and y values with the colormap set with `cmap`.
`fig, ax = plt.subplots(nrows, ncols)`	Tells python to create a grid of subplots with `nrows` number of rows and `ncols` number of columns. `ax` can be indexed to refer to a specific subplot, f.ex. the subplot on the first row and second column would be `ax[0,1]`.
`gaussian_filter1d(data, sigma = filter_width)`	Apply a 1D gaussian filter that smoothes the data to reduce the prominence of noise in a plot. `sigma = filter_width` designates the width of the filter.
`ax_twin = ax1.twinx()`	Makes an extra axis in a subplot (the subplot of `ax1`) such that data with different units can be plotted together in the same plot.

Exercise: Run the cell below to get the CSD and the population firing rates.

# provided

lfp_neo = utils.create_neo_lfp(lfp_values=lfp_trial_avg.values.T, sampling_frequency=sampling_frequency_lfp) # The lfp must be transposed to have time on the first dimension for neo
coords = utils.add_um_units(channels_depth_vals=channels_depth.values)
csd_neo = estimate_csd(lfp_neo, coordinates=coords, method='DeltaiCSD', diam=400*pq.um).T # The csd must be transposed to get channels on first dimension for the plot

pop_firing_rate_V1_trial_avg = np.mean(spike_counts_V1, axis=(0,1))*1E3 # multiply with 1E3 to scale to Hz

# get sampling time points of LFP and spike data
time_lfp = lfp_trial_avg.time_lfp
time_spikes = spike_counts_V1.time_spikes

discrete filter coefficients: 
b = [ 0.607 1.000 0.607 ],                
a = [ 2.213 ]

Example: Plot the firing rate of V1 and the CSD from V1 together in the same plot using ax.twinx().

fig, ax = plt.subplots(nrows = 1, figsize = (4,4))

ax1 = ax
im = ax1.pcolormesh(time_lfp, channels_depth, csd_neo, cmap = 'RdBu_r')
ax1.set(xlabel='Time from stimulus onset (ms)', ylabel=r'Depth ($\mu$m)')

ax2 = ax1.twinx()
ax2.plot(time_spikes, pop_firing_rate_V1_trial_avg)
ax2.set(ylabel='Firing rate (Hz)');

Exercise: Plot the population firing rate of LGN and the CSD from V1 together in the same plot using ax.twinx(). You can apply a gaussian filter to the firing rate if you want to reduce the noise in the plot.

# provided

pop_firing_rate_LGN_trial_avg = spike_counts_LGN.mean(['trial_nr', 'unit_id_LGN'])*1E3 # multiply with 1E3 to scale to Hz

Solution

# solution
fig, ax = plt.subplots(nrows = 1, figsize = (4,4))

ax1 = ax
im = ax1.pcolormesh(time_lfp, channels_depth, csd_neo, cmap = 'RdBu_r')
ax1.set(xlabel='Time from stimulus onset (ms)', ylabel=r'Depth ($\mu$m)')

ax2 = ax1.twinx()
ax2.plot(time_spikes, pop_firing_rate_LGN_trial_avg)
ax2.set(ylabel='Firing rate (Hz)');

Exercise: Plot the population firing rate of LGN and the CSD from V1 for the first 100 ms after stimulus onset using ax.twinx().

# provided

time_window_start = 0
time_window_end = 100

# get csd and time points of lfp samples in time window
mask_window_lfp = (time_lfp > time_window_start) & (time_lfp < time_window_end)
time_window_lfp = time_lfp[mask_window_lfp]
csd_window = csd_neo[:, mask_window_lfp]

# get population firing rate and time points of spike data samples in time window
mask_window_spikes = (time_spikes > time_window_start) & (time_spikes < time_window_end)
time_window_spikes = time_spikes[mask_window_spikes]
pop_firing_rate_LGN_window = pop_firing_rate_LGN_trial_avg[mask_window_spikes]

Solution

# solution
fig, ax = plt.subplots(nrows = 1, figsize = (4,4))

ax1 = ax
im = ax1.pcolormesh(time_lfp, channels_depth, csd_neo, cmap = 'RdBu_r')
ax1.set(xlabel='Time from stimulus onset (ms)', ylabel=r'Depth ($\mu$m)')

ax2 = ax1.twinx()
ax2.plot(time_spikes, pop_firing_rate_LGN_trial_avg)
ax2.set(ylabel='Firing rate (Hz)');

Exercise: Plot the firing rate of all three structures (LGN, V1, and LM) and the CSD from V1 together in the same plot using ax.twinx(). Normalize the firing rates with the max value. If you want, you can smooth the firing rates with a gaussian filter.

# provided

pop_firing_rate_LM_trial_avg = np.mean(spike_counts_LM, axis=(0,1))*1E3 # multiply with 1E3 to scale to Hz
pop_firing_rate_LM_window = pop_firing_rate_LM_trial_avg[mask_window_spikes]

pop_firing_rate_V1_trial_avg = np.mean(spike_counts_V1, axis=(0,1))*1E3
pop_firing_rate_V1_window = pop_firing_rate_V1_trial_avg[mask_window_spikes]

Solution

fig, ax = plt.subplots(nrows = 1, figsize = (4,4))

ax1 = ax
im = ax1.pcolormesh(time_window_lfp, channels_depth, csd_window, cmap = 'RdBu_r')
ax1.set(xlabel='Time from stimulus onset (ms)', ylabel=r'Depth ($\mu$m)')

ax2 = ax1.twinx()
ax2.plot(time_window_spikes, pop_firing_rate_LGN_window)
ax2.plot(time_window_spikes, pop_firing_rate_V1_window)
ax2.plot(time_window_spikes, pop_firing_rate_LM_window)

ax2.set(ylabel='Firing rate (Hz)')

Section 4: Interpolating Signals with Different Sampling Rates to Quantify Their Association

Since the spikes and the LFP are sampled at different frequencies (here, the spikes have a sampling frequency of 1000 Hz and the LFP has a sampling frequency of 1250 Hz), we have to interpolate the signals onto the same basis to be able to calculate statistical associations between them.

Code	Description
`f = RectBivariateSpline(x,y,z)`	Create function to use for interpolation of signal.
`signal_interpolated = f(x_new, y_new)`	Create the interpolated signal using the function created with the code above and the new x- and y-basis as input.
`correlation, p_value = scipy.stats.pearsonr(signalA, signalB)`	Calculate the pearson correlation between signal A and signal B.

Exercises

from scipy.interpolate import RectBivariateSpline

from scipy.stats import pearsonr

Run cell below to calculate the CSD and trial-averaged population firing rate in V1.

# provided

lfp_neo = utils.create_neo_lfp(lfp_values=lfp_trial_avg.values.T, sampling_frequency=sampling_frequency_lfp) # The lfp must be transposed to have time on the first dimension for neo
coords = utils.add_um_units(channels_depth_vals=channels_depth.values)
csd_neo = estimate_csd(lfp_neo, coordinates=coords, method='DeltaiCSD', diam=400*pq.um).T # The csd must be transposed to get channels on first dimension for the plot

# calculate population firing rate in V1
pop_firing_rate_V1_trial_avg = spike_counts_V1.mean(['unit_id_V1', 'trial_nr'])*1E3 # multiply with 1E3 to scale to Hz

# get sampling time points of LFP and spike data
time_lfp = lfp_trial_avg.time_lfp
time_spikes = spike_counts_V1.time_spikes
csd_neo

discrete filter coefficients: 
b = [ 0.607 1.000 0.607 ],                
a = [ 2.213 ]

AnalogSignal with 1875 channels of length 23; units A/m**2; datatype float64
annotations: {'x_coords': array([[  0.],
         [ 40.],
         [ 80.],
         [120.],
         [160.],
         [200.],
         [240.],
         [280.],
         [320.],
         [360.],
         [400.],
         [440.],
         [480.],
         [520.],
         [560.],
         [600.],
         [640.],
         [680.],
         [720.],
         [760.],
         [800.],
         [840.],
         [880.]]) * um}
sampling rate: 1250.0 Hz
time: 0.0 s to 0.0184 s

Example: Interpolate the CSD for the first 100 ms after flash onset so that it has the same basis (sampling rate) as the firing rate in LGN Make a figure with subplots where you plot both the original and the interpolated CSD. The CSD should be have the same main pattern, but the ticks on the x-axis should be different.

time_basis_csd = np.linspace(0., 1., csd_neo.shape[1])
channel_basis_csd = np.linspace(0., 1., csd_neo.shape[0])
interpolation_function = RectBivariateSpline(time_basis_csd, channel_basis_csd, csd_neo.T)

new_time_basis = np.linspace(0., 1., len(pop_firing_rate_LGN_trial_avg))
new_channel_basis = channel_basis_csd # only interpolate on time axis, so channel basis is unchanged
csd_interpolated = interpolation_function(new_time_basis, new_channel_basis).T

# do a visualization check to see that the CSD is still the same after interpolation

fig, ax = plt.subplots(nrows = 2, figsize=(3,6))
ax[0].imshow(csd_neo, cmap='RdBu_r', aspect='auto')
ax[0].set_title('CSD in V1 before interpolation')

ax[1].imshow(csd_interpolated, cmap='RdBu_r', aspect='auto')
ax[1].set_title('CSD in V1 after interpolation')

fig.tight_layout()

Exercise: Select the time window from 0 to 50 milliseconds after stimulus onset and interpolate the CSD so that the CSD and the firing rate in LGN have the same time basis (the same sampling rate). Plot original and interpolated CSD in this 250 ms time window in different subplots to do a visual check of the interpolation.

# provided

time_window_start = 0
time_window_end = 50

# get csd and time points of lfp samples in time window
mask_window_lfp = (time_lfp > time_window_start) & (time_lfp < time_window_end)
time_window_lfp = time_lfp[mask_window_lfp]
csd_window = csd_neo[:, mask_window_lfp]

# get population firing rate and time points of spike data samples in time window
mask_window_spikes = (time_spikes > time_window_start) & (time_spikes < time_window_end)
time_window_spikes = time_spikes[mask_window_spikes]
pop_firing_rate_LGN_window = pop_firing_rate_LGN_trial_avg[mask_window_spikes]

Solution

time_basis_csd = np.linspace(0., 1., csd_window.shape[1])
channel_basis_csd = np.linspace(0., 1., csd_window.shape[0])

new_time_basis = np.linspace(0., 1., len(pop_firing_rate_LGN_window))
new_channel_basis = channel_basis_csd # only interpolate on time axis, so channel basis is unchanged

interpolation_function = RectBivariateSpline(time_basis_csd, channel_basis_csd, csd_window.T)

csd_interpolated = interpolation_function(new_time_basis, new_channel_basis).T

Solution

# do a visualization check to see that the CSD is still the same after interpolation

fig, ax = plt.subplots(nrows = 2, figsize=(3,6))
ax[0].imshow(csd_window, cmap='RdBu_r', aspect='auto')
ax[0].set_title('CSD in V1 before interpolation')

ax[1].imshow(csd_interpolated, cmap='RdBu_r', aspect='auto')
ax[1].set_title('CSD in V1 after interpolation')

fig.tight_layout()

When the firing rate and the CSD have the same basis (an equal number of samples in a time window), it’s possible to calculate the statistical association between them.

Exercise: Calculate the correlation between the interpolated CSD from channel index 8 and the firing rate in LGN in the 50 ms after stimulus onset time window. (Remember that since current sinks have negative values, the sign of the correlation would be negative if a peak in the LGN firing rate and a current sink coincide in time.)

Solution

ichan = 8
pearsonr(csd_interpolated[ichan], pop_firing_rate_LGN_window)

PearsonRResult(statistic=-0.4742870888496559, pvalue=0.000502171851041129)

Exercise: Plot the CSD trace in channel 8 and the population firing rate in LGN together in the same plot using ax.twinx().

Solution

fig, ax = plt.subplots()

ax1 = ax
ax1.plot(pop_firing_rate_LGN_window, label = 'Firing rate LGN', color = 'r')

ax2=ax1.twinx()
ax2.plot(csd_interpolated[ichan], label = f'CSD channel {ichan}')
fig.legend();

ax1.set_xlabel('Time from stimulus onset (ms)');