CSD Analysis Theory and Methods

The LFP can be challenging to interpret in terms of the underlying neurophysiology. One method that helps make the LFP more interpretable is current source density (CSD) analysis. The CSD is calculated from the LFP and shows where the current sinks and sources generating the LFP are located in space and when they arise. A current sink shows where positive ions enter a neuron/population of neurons, while a current source shows where positive ions leave a neuron/population of neurons.

When a neuron receives excitatory input, positive ions flush into the neuron at the synapse. This will show up as a current sink in a CSD plot. Therefore, a current sink can tell you where and when a population of neurons gets input. There are multiple different neuronal processes that can give rise to a current sink or a current source, so the interpretation is still not always straightforward, but the CSD at least brings us one step closer to interpreting the electrical signals we record in terms of the neurophysiology.

Setup

Import Libraries

import numpy as np
import neo
import quantities as pq
from scipy.stats import pearsonr
import matplotlib.pyplot as plt
from matplotlib.colors import TwoSlopeNorm
from elephant.current_source_density import estimate_csd

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.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_data_time_window(data, time, time_window_start, time_window_end, is_trial_avg = True):
    """Get lfp or csd data in selected time window"""
    mask_time_window = (time > time_window_start) & (time < time_window_end)

    # get data in time window
    if is_trial_avg:
        data_window = data[:, mask_time_window]
    else:
        data_window = data[:, :, mask_time_window]

    return data_window

def correlate_pairwise_trials(csd):
    '''Correlate csd from all trials pairwise'''
    corr_csd_all_trials = []

    for itrialA in range(csd.shape[1]):
        corr_csd_temp = []
        for itrialB in range(csd.shape[1]):
            corr_csd_temp.append(pearsonr(csd[:,itrialA].flatten(), csd[:,itrialB].flatten())[0])
        corr_csd_all_trials.append(corr_csd_temp)

    return np.array(corr_csd_all_trials)

def euclidian_norm_pairwise_trials(csd):
    '''Correlate csd from all trials pairwise'''
    eucnorm_csd_all_trials = []

    for itrialA in range(csd.shape[1]):
        eucnorm_csd_temp = []
        for itrialB in range(csd.shape[1]):
            eucnorm_csd_temp.append(np.linalg.norm(csd[:,itrialA]-csd[:,itrialB]))
        eucnorm_csd_all_trials.append(eucnorm_csd_temp)

    return np.array(eucnorm_csd_all_trials)


def plot_csd(csd):
    plt.figure()
    plt.imshow(
        np.array(csd.T),                          # The csd must be transposed to get channels on first dimension for the plot
        aspect = 'auto', 
        cmap = 'RdBu_r',  # 
        norm=TwoSlopeNorm(vcenter=0),
    )  # Puts white a the zero value, so blue is always negative and red is always positive
    plt.ylabel('Channel nr.')
    plt.colorbar();


def calc_csd(lfp, sampling_freq, channel_depths, method='StandardCSD', **csd_kwargs) -> 'neo':
    lfp_neo = create_neo_lfp(lfp, sampling_frequency=sampling_freq) # The lfp must be transposed to have time on the first dimension for neo
    coords = add_um_units(channels_depth_vals=channel_depths.values)  # create array with channel positions and units
    csd_neo = estimate_csd(lfp_neo, coordinates=coords, method=method, **csd_kwargs)
    return csd_neo

class utils:
    correlate_pairwise_trials = correlate_pairwise_trials
    euclidian_norm_pairwise_trials = euclidian_norm_pairwise_trials
    calc_csd = calc_csd
    plot_csd = plot_csd

Download dataset

import owncloud

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

True

Section 1: Calculate the Current Source Density (CSD) from the Recorded LFP Using the “Standard” Method.

In this section, you will learn how to calculate the CSD from the LFP using the so-called standard method.

Load Dataset

Run the cell below to load dataset and extract LFP data.

Exercises

import xarray as xr
from pathlib import Path

dataset_path = Path('dataset_session_754312389.nc')

dataset = xr.load_dataset(dataset_path)
dataset = dataset[[
    'lfp_V1',
    'csd_V1'
]]
dataset

<xarray.Dataset> Size: 52MB
Dimensions:        (channel_depth: 23, trial_nr: 75, time_lfp: 1875,
time_csd: 1875)
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


time_csd       (time_csd) float64 15kB -1e+03 -999.2 -998.4 … 498.4 499.2
Data variables:
lfp_V1         (channel_depth, trial_nr, time_lfp) float64 26MB -8.516e-0…
csd_V1         (channel_depth, trial_nr, time_csd) float64 26MB -0.04729 …
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

Example: Calculate and plot the CSD for the first trial (trial index 0). You can get the lfp from the first trial using traditional array indexing (e.g. lfp_V1[:, 0]) or the sel method (e.g. lfp_V1.sel(trial_nr = 0)).

lfp = dataset['lfp_V1'].sel(trial_nr = 0)
csd = utils.calc_csd(lfp, dataset.sampling_frequency_lfp, dataset.channel_depth, method='StandardCSD')
utils.plot_csd(csd)

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

Exercise: Calculate and plot the CSD for trial at index 30.

Solution

lfp = dataset['lfp_V1'].sel(trial_nr = 30)
csd = utils.calc_csd(lfp, dataset.sampling_frequency_lfp, dataset.channel_depth, method='StandardCSD')
utils.plot_csd(csd)

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

Exercise: Calculate and plot the CSD for the trial averaged LFP (Tip: lfp.mean(dim='trial_nr'))

Solution

lfp = dataset['lfp_V1'].mean(dim = 'trial_nr')
csd = utils.calc_csd(lfp, dataset.sampling_frequency_lfp, dataset.channel_depth, method='StandardCSD')
utils.plot_csd(csd)

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

Example: Using the gaussian_filter() function, smooth the CSD from across depth (over the channels) with a gaussian filter with sigma = (0,1). (The 0 means that no smoothing is across over the time dimension.) Plot the CSD after smoothing with the pcolormesh function.

from scipy.ndimage import gaussian_filter

lfp = dataset['lfp_V1'].mean(dim='trial_nr')
csd = utils.calc_csd(lfp, dataset.sampling_frequency_lfp, dataset.channel_depth)
csd_smoothed = gaussian_filter(csd, sigma=(1, 0))
utils.plot_csd(csd_smoothed)

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

Exercise: Smooth the CSD across depth (over the channels) with a gaussian filter with sigma = (0,4). Plot the CSD after smoothing.

What happened to the minor and the major sinks and sources in the plot when you compare to CSD plots with no or less smoothing?

Solution

lfp = dataset['lfp_V1'].mean(dim='trial_nr')
csd = utils.calc_csd(lfp, dataset.sampling_frequency_lfp, dataset.channel_depth)
csd_smoothed = gaussian_filter(csd, sigma=(0, 4))
utils.plot_csd(csd_smoothed)

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

Exercise: Smooth the CSD across time with a gaussian filter with sigma = (4,0). Plot the CSD after smoothing.

Solution

lfp = dataset['lfp_V1'].mean(dim='trial_nr')
csd = utils.calc_csd(lfp, dataset.sampling_frequency_lfp, dataset.channel_depth)
csd_smoothed = gaussian_filter(csd, sigma=(4, 0))
utils.plot_csd(csd_smoothed)

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

Exercise: Smooth the CSD across depth and time with a gaussian filter with sigma = (2,2). Plot the CSD after smoothing.

Solution

lfp = dataset['lfp_V1'].mean(dim='trial_nr')
csd = utils.calc_csd(lfp, dataset.sampling_frequency_lfp, dataset.channel_depth)
csd_smoothed = gaussian_filter(csd, sigma=(2, 2))
utils.plot_csd(csd_smoothed)

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

Comment: While smoothing and interpolation can improve your visalization by making your plots more readable, excessive smoothing can make them misleading, so use with care.

Section 2: Slicing Time Points, Plotting Time, and Making 2D Plots Using plt.pcolormesh()

Exercises

Example: Extract and plot the CSD calculated from the trial averaged LFP from stimulus onset to the end of the trial.

mask = dataset.time_lfp > 0
lfp = dataset['lfp_V1'].where(mask, drop=True).mean(dim='trial_nr')
csd = utils.calc_csd(lfp, dataset.sampling_frequency_lfp, dataset.channel_depth, method='StandardCSD')
utils.plot_csd(csd)

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

Exercise: Extract and plot the CSD calculated from the trial averaged LFP from the beginning of the trial to the stimulus onset.

Solution

mask = dataset.time_lfp > -1000
lfp = dataset['lfp_V1'].where(mask, drop=True).mean(dim='trial_nr')
csd = utils.calc_csd(lfp, dataset.sampling_frequency_lfp, dataset.channel_depth, method='StandardCSD')
utils.plot_csd(csd)

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

Exercise: Extract and plot the CSD from stimulus onset until 100 msec after stimulus onset.

Approximately when does the first current sink (blue blob) arise? Is the position of this blob between channels 5 and 10 - roughly layer 4 of the primary visual cortex? If so, the CSD plot is correct.

Useful Code:

mask = (0 < dataset['time_lfp']) & (dataset['time_lfp'] < 100)

Solution

mask = (0 < dataset.time_lfp) & (dataset.time_lfp < 100)
lfp = dataset['lfp_V1'].where(mask, drop=True).mean(dim='trial_nr')
csd = utils.calc_csd(lfp, dataset.sampling_frequency_lfp, dataset.channel_depth)
utils.plot_csd(csd)

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

Example: Plot the trial averaged CSD in the first 100 ms after stimulus onset using plt.pcolormesh function. Pass the time points that the LFP was sampled as the first parameter and the depth of channels (channels_depth) as the second parameter to the pcolormesh function.

mask = (0 < dataset.time_lfp) & (dataset.time_lfp < 100)
lfp = dataset['lfp_V1'].where(mask, drop=True).mean(dim='trial_nr')
csd = utils.calc_csd(lfp, dataset.sampling_frequency_lfp, dataset.channel_depth)

plt.figure()
plt.pcolormesh(lfp.time_lfp, dataset.channel_depth, csd.T, cmap = 'RdBu_r')
plt.xlabel('Time (ms)')
plt.ylabel(r'Depth ($\mu$m)');

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

Exercise: Plot the trial averaged CSD in the first 250 ms after stimulus onset using plt.pcolormesh function. Pass the time points that the LFP was sampled as the first parameter and the depth of channels (channels_depth) as the second parameter to the pcolormesh function.

Solution

mask = (0 < dataset.time_lfp) & (dataset.time_lfp < 250)
lfp = dataset['lfp_V1'].where(mask, drop=True).mean(dim='trial_nr')
csd = utils.calc_csd(lfp, dataset.sampling_frequency_lfp, dataset.channel_depth)

plt.figure()
plt.pcolormesh(lfp.time_lfp, dataset.channel_depth, csd.T, cmap = 'RdBu_r')
plt.xlabel('Time (ms)')
plt.ylabel(r'Depth ($\mu$m)');

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

Exercise: Get the CSD from the from 200 ms before stimulus onset to 200 ms after stimulus onset and plot it using plt.pcolormesh function.

Solution

mask = (-200 < dataset.time_lfp) & (dataset.time_lfp < 200)
lfp = dataset['lfp_V1'].where(mask, drop=True).mean(dim='trial_nr')
csd = utils.calc_csd(lfp, dataset.sampling_frequency_lfp, dataset.channel_depth)

plt.figure()
plt.pcolormesh(lfp.time_lfp, dataset.channel_depth, csd.T, cmap = 'RdBu_r')
plt.xlabel('Time (ms)')
plt.ylabel(r'Depth ($\mu$m)');

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

Section 3: Calculate the Current Source Density (CSD) from the Recorded LFP Using the Delta Source Method.

The method you used above to compute the CSD is called the “standard” method simply because it was developed first. The standard method assumes that the CSD is infinitely constant lateral to the probe, which does not always hold true. Sometimes, this approximation is good enough to not cause any major errors in the CSD estimate, but in other cases it can result in an incorrect CSD amplitude and to spurious current sinks and sources in the estimate.

One alternative method that has been developed to deal with this issue is called the “delta source” method, which only assumes that the CSD is constant within a certain diameter lateral to the probe. In this section, you will explore the calculation of CSD using the delta source method.

Exercises

Example: Calculate and plot the CSD for the trial averaged LFP using Delta iCSD method rather than the standard method.

Set the method parameter in the utils.calc_csd function to DeltaiCSD. Set the diam parameter to 800*pq.um.

diameter = 800*pq.um
lfp = dataset['lfp_V1'].mean(dim='trial_nr')
csd = utils.calc_csd(lfp, dataset.sampling_frequency_lfp, dataset.channel_depth, method='DeltaiCSD', diam=diameter)
utils.plot_csd(csd)

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

Exercise: Change the diameter parameter to 200 $\mu m$ and calculate and plot the the CSD again. What changed?

Solution

diameter = 200*pq.um
lfp = dataset['lfp_V1'].mean(dim='trial_nr')
csd = utils.calc_csd(lfp, dataset.sampling_frequency_lfp, dataset.channel_depth, method='DeltaiCSD', diam=diameter)
utils.plot_csd(csd)

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

Exercise: Change the diameter parameter to 1600 $\mu m$ and calculate and plot the the CSD again. What changed?

Solution

diameter = 1600*pq.um
lfp = dataset['lfp_V1'].mean(dim='trial_nr')
csd = utils.calc_csd(lfp, dataset.sampling_frequency_lfp, dataset.channel_depth, method='DeltaiCSD', diam=diameter)
utils.plot_csd(csd)

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

Exercise: Calculate the CSD with the standard method (method='StandardCSD') and plot it again.

Which of the three plots above (where the CSD was calculated with the delta method with different diameters) does this CSD plot resemble the most? Is it the plot you would expect, considering the diameter used with the delta method in that plot?

Solution

lfp = dataset['lfp_V1'].mean(dim='trial_nr')
csd = utils.calc_csd(lfp, dataset.sampling_frequency_lfp, dataset.channel_depth, method='StandardCSD')
utils.plot_csd(csd)

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

Comment: As you saw above, the diameter parameter you set for the delta method will affect the amplitude and the distribution of sinks and sources in the CSD estimate. Usually, one won’t exactly know the correct value for diameter. You have to settle for a reasonable estimate based on the physiology of the structure in question. For the cortex, a reasonable value to use is the diameter of a cortical column, which is approximately 300-500 micrometer in mice.

Section 4: Metrics to Quantify CSD Similarity: Correlation, Euclidian Distance, and Wasserstein Distance.

Sometimes you want to compare the CSD from different trials or different animals. Visualization is a great way to do that, but you may need to quantify the level of similarity between CSD plots too. Below are a few methods that can be used for that.

Exercises

from scipy.stats import wasserstein_distance, pearsonr

Run the cell below to get pre-calculated CSD from all trials and select the CSD from the first 100 ms after stimulus onset.

time_window_start = 0 # at stimulus onset
time_window_end = 100

# load precomputed csd in all trials form dataset
csd = dataset['csd_V1']
time_csd = csd.time_csd
mask_time_window = np.logical_and(time_window_start < time_csd, 
                                  time_window_end > time_csd)

csd = csd[:,:,mask_time_window]
csd.shape

(23, 75, 125)

Example: Make a scatter plot of the csd values in trial 20 and trial 21

# provided
csd_trialA = csd.sel(trial_nr = 20).values.flatten()
csd_trialB = csd.sel(trial_nr = 21).values.flatten()

plt.scatter(csd_trialA, csd_trialB, marker='.', s=1.5);

Exercise: Calculate and display the correlation between the CSD from trial 20 and trial 21.

pearsonr(csd_trialA, csd_trialB)

PearsonRResult(statistic=0.6614942496566328, pvalue=0.0)

Exercise: Calculate the correlation between the CSD from trial 0 and the last trial (set itrialB to 74).

Are these trials more or less similar than the neighboring trials in the previous exercise? Make a scatter plot where you plot the values of the CSD in trial 0 against the values of the CSD in trial 74 to check that the result makes sense.

Solution

csd_trialA = csd.sel(trial_nr = 0).values.flatten()
csd_trialB = csd.sel(trial_nr = 74).values.flatten()

plt.scatter(csd_trialA, csd_trialB, marker='.', s=1.5);
pearsonr(csd_trialA, csd_trialB)

PearsonRResult(statistic=0.10428248676854059, pvalue=2.0900224966659755e-08)

Exercise: Calculate the correlation between the CSD from the first trial and the trial averaged CSD.

pearsonr(csd_trialA.flatten(), np.mean(csd, axis=1).values.flatten())

PearsonRResult(statistic=0.24234596134818845, pvalue=1.0600984623522463e-39)

Demo: Run the cell below to calculate and plot the pairwise correlation between the CSD from all trials.

Does this plot suggest that neighboring trials are typically more correlated?

corr_csd_all_trials = utils.correlate_pairwise_trials(csd.values)

plt.imshow(corr_csd_all_trials, cmap = 'bwr')
plt.xlabel('Trial nr.')
plt.ylabel('Trial nr.')
plt.colorbar();

Exercise: Calculate the euclidian distance between the CSD from the first and the second trial.

csd_trialA = csd.sel(trial_nr = 0).values.flatten()
csd_trialB = csd.sel(trial_nr = 1).values.flatten()

np.linalg.norm(csd_trialA - csd_trialB)

6.763006171185026

Exercise: Calculate the euclidian distance between the CSD from trial 0 and the last trial.

csd_trialA = csd.sel(trial_nr = 0).values.flatten()
csd_trialB = csd.sel(trial_nr = 74).values.flatten()

np.linalg.norm(csd_trialA - csd_trialB)

7.175334725036059

Exercise: Calculate and plot the pairwise euclidian distance between the the CSD from all trials.

An important difference between the euclidian distance and the correlation is that amplitude differences also matters for the euclidian distance, so this plot will not necessarily have the same structure as the correlation plot above.

Solution

csd_pairwise_euclid = utils.euclidian_norm_pairwise_trials(csd)

plt.imshow(csd_pairwise_euclid, cmap = 'bwr')
plt.xlabel('Trial nr.')
plt.ylabel('Trial nr.')
plt.colorbar();

Exercise: Calculate the Wasserstein distance between the CSD from the first and the second trial.

csd_trialA = csd.sel(trial_nr = 0).values.flatten()
csd_trialB = csd.sel(trial_nr = 1).values.flatten()

wasserstein_distance(csd_trialA.flatten(), csd_trialB.flatten())

0.04407472847070443

Exercise: Calculate the Wasserstein distance between the CSD from the first and the last trial.

csd_trialA = csd.sel(trial_nr = 0).values.flatten()
csd_trialB = csd.sel(trial_nr = 74).values.flatten()

wasserstein_distance(csd_trialA.flatten(), csd_trialB.flatten())

0.03141649491519594

Discussion Exercise:

• What are some potential reasons for using one or another metric to quantify similarity between two CSD plots?