LFP Analysis

LFP Analysis#

Tutorial overview#

This Jupyter notebook will demonstrate how to access and analyze LFP data from the Neuropixels Visual Behavior dataset. LFP, which stands for “local field potential,” contains information about low-frequency (0.1-500 Hz) voltage fluctations around each recording site. It’s complementary to the spiking activity, and can be analyzed on its own or in conjunction with spikes.

This tutorial will cover the following topics:

Selecting probes to analyze
Loading LFP data
Aligning data in time
Aligning data in space
Exploring pre-computed CSD plots

This tutorial assumes you’ve already created a data cache, or are working with the files on AWS. If you haven’t reached that step yet, we recommend going through the data access tutorial first.

Functions related to analyzing spike data will be covered in other tutorials. For a full list of available tutorials, see the SDK documentation.

Learning about the available probes#

First let’s import the VisualBehaviorNeuropixelsProjectCache and a few other helpful packages:

from allensdk.brain_observatory.behavior.behavior_project_cache.\
    behavior_neuropixels_project_cache \
    import VisualBehaviorNeuropixelsProjectCache

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

/opt/envs/allensdk/lib/python3.8/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

cache_dir = "/data/"

Now let’s instantiate a cache from the Amazon S3 bucket. The cache_dir should be a path specifying where to download data on your local system:

cache = VisualBehaviorNeuropixelsProjectCache.from_local_cache(
            cache_dir=cache_dir, use_static_cache=True)

Finding LFP data of interest

We can look at the cache probes table to identify which probe insertions have valid LFP data

probes = cache.get_probe_table()
valid_lfp = probes[probes['has_lfp_data']]

print('Fraction of insertions with valid LFP: ', len(valid_lfp)/len(probes))

Fraction of insertions with valid LFP:  0.9602209944751381

Now let’s find a probe insertion with good LFP data that passed through VISp.

valid_lfp[valid_lfp['structure_acronyms'].str.contains("'VISp',")].head()

	ecephys_session_id	name	sampling_rate	lfp_sampling_rate	phase	has_lfp_data	unit_count	channel_count	structure_acronyms	file_id
ecephys_probe_id
1044506934	1044385384	probeC	30000.049852	2500.004154	1.0	True	307	384	['MB', 'MRN', 'POST', 'SCig', 'VISp', 'root']	1
1044520064	1044389060	probeC	30000.146997	2500.012250	1.0	True	417	384	['MB', 'NB', 'POL', 'POST', 'PRE', 'SUB', 'TH'...	7
1044791095	1044594870	probeC	30000.054069	2500.004506	1.0	True	365	384	['HPF', 'MB', 'MRN', 'POST', 'SCig', 'VISp', '...	12
1044827086	1044597824	probeC	30000.148958	2500.012413	1.0	True	378	384	['DG', 'MB', 'PIL', 'POST', 'PRE', 'SGN', 'SUB...	18
1046469927	1046166369	probeC	30000.150477	2500.012540	1.0	True	519	384	['HPF', 'MB', 'MRN', 'POST', 'SCig', 'VISp', '...	24

Below we’ll select an insertion from this table (1064735073) and grab the session data for it. We can get the session ID from the probes table:

session_id = probes.loc[1064735073]['ecephys_session_id']
session = cache.get_ecephys_session(
            ecephys_session_id=session_id)

Aligning LFP data to a stimulus#

In the above example, we selected LFP data based on an arbitrary time span (100 to 101 seconds). For many analyses, however, you’ll want to align the data to the onset of a particular type of stimulus.

First, we need to select some stimulus presentations to use. Below, we’ll use the full field flashes presented in stimulus block 4 of the Visual Behavior Neuropixels experiments. To learn more about the stimuli used in these experiments, check out the documentation

stim_presentations = session.stimulus_presentations
flashes = stim_presentations[stim_presentations['stimulus_name'].str.contains('flash')]
presentation_times = flashes.start_time.values
presentation_ids = flashes.index.values

First, let’s make a convenience function that helps us align the LFP to times of interest. Because we’re using xarrays, the alignment operation is fast, and doesn’t require any for loops! There’s a lot going on here, so we recommend referring to the pandas and xarray documentation if anything is confusing:

def align_lfp(lfp, trial_window, alignment_times, trial_ids = None):
    '''
    Aligns the LFP data array to experiment times of interest
    INPUTS:
        lfp: data array containing LFP data for one probe insertion
        trial_window: vector specifying the time points to excise around each alignment time
        alignment_times: experiment times around which to excise data
        trial_ids: indices in the session stim table specifying which stimuli to use for alignment.
                    None if aligning to non-stimulus times
    
    OUTPUT:
        aligned data array with dimensions channels x trials x time
    '''
    
    time_selection = np.concatenate([trial_window + t for t in alignment_times])
    
    if trial_ids is None:
        trial_ids = np.arange(len(alignment_times))
        
    inds = pd.MultiIndex.from_product((trial_ids, trial_window), 
                                      names=('presentation_id', 'time_from_presentation_onset'))

    ds = lfp.sel(time = time_selection, method='nearest').to_dataset(name = 'aligned_lfp')
    ds = ds.assign(time=inds).unstack('time')

    return ds['aligned_lfp']

aligned_lfp = align_lfp(lfp, np.arange(-0.5, 0.5, 1/500), presentation_times, presentation_ids)

aligned_lfp is a DataArray with dimensions of channels x trials x time. It’s been downsampled to 500 Hz by changing the time step in the trial_window argument of the align_lfp function.

Note that we can get the channels IDs for each channel in this DataArray. Let’s use the session channels table to map these to the probe and mark the surface of the brain.

chans = session.get_channels()
lfp_chan_depths = [chans.loc[c]['probe_vertical_position'] for c in lfp.channel.values]

chans_in_brain = chans[(chans['probe_id']==1064735073)&(~chans['structure_acronym'].str.contains('root'))]
first_channel_in_brain_position = chans_in_brain['probe_vertical_position'].max()

fig, ax = plt.subplots()
fig.suptitle('Flash aligned mean LFP')
im = ax.pcolor(aligned_lfp.time_from_presentation_onset.values, lfp_chan_depths, aligned_lfp.mean(dim='presentation_id').data)
_ = plt.colorbar(im, fraction=0.036, pad=0.04)
_ = plt.xlabel('Time from flash onset (s)')
_ = plt.ylabel('Channel Position from Tip (um)')

ax.axvline(0, c='w', ls='dotted')
ax.axvline(0.25, c='w', ls='dotted')
ax.axhline(first_channel_in_brain_position, c='w')
ax.text(-0.4, first_channel_in_brain_position+50, 'brain surface', c='w')

/tmp/ipykernel_5031/759061444.py:3: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.
  im = ax.pcolor(aligned_lfp.time_from_presentation_onset.values, lfp_chan_depths, aligned_lfp.mean(dim='presentation_id').data)

Text(-0.4, 3570, 'brain surface')

../_images/69bdea3c340edabb3763c9778b646e2fe5cc7535c8f53b966c1038df1cb23fb6.png

Here we see the effect of a 250 ms flash stimulus on the LFP. There are two large responses in cortex (the first ~700 um below the brain surface), one corresponding to the stimulus onset (around 50 ms), and one corresponding to the stimulus offset (around sample 300 ms).

You can use the code sample above to align the LFP to any type of event (e.g. spike times, running onset, optogenetic stimuli) just by changing the trial_window and time_selection variables.

Aligning LFP data to units#

The previous section demonstrated how to align the LFP in time. What if we want to extract the LFP at a particular location in space, corresponding to the location of a unit we’re analyzing?

Let’s start by finding a well-isolated unit whose peak channel is included in our LFP data.

Once we’ve selected a unit of interest, we can align the LFP data to its spike times:

sess_units  = session.get_units()

#Grab units whose peak channels are in the LFP data, have relatively low isi violations and high amplitude spikes
units_on_lfp_chans = sess_units[(sess_units.peak_channel_id.isin(lfp.channel.values)) &
                                (sess_units.isi_violations < 0.5) &
                                (sess_units.amplitude > 200)]

#Merge this curated unit table with the channel table to get CCF locations for these units
units_on_lfp_chans = units_on_lfp_chans.merge(chans, left_on='peak_channel_id', right_index=True)

#Select a unit in V1
v1_units = units_on_lfp_chans[units_on_lfp_chans.structure_acronym.str.contains('VISp')]
unit_id = v1_units.index.values[5]

#Get the peak channel ID for this unit (the channel on which it had the greatest spike amplitude)
peak_chan_id = units_on_lfp_chans.loc[unit_id]['peak_channel_id']
peak_probe_position = units_on_lfp_chans.loc[unit_id]['probe_vertical_position']

Using unit_id and peak_chan_id, we can select the spikes and LFP within an arbitrary time interval. Note that we can use method='nearest' when selecting the LFP data channel if our peak channel isn’t in the LFP DataArray. Here this is unnecessary since we’ve filtered for units with peak channels in the LFP data.

start_time = 500
end_time = 510

spike_times = session.spike_times[unit_id]

times_in_range = spike_times[(spike_times > start_time) & (spike_times < end_time)]

lfp_data = lfp.sel(time = slice(start_time, end_time))
lfp_data = lfp_data.sel(channel = peak_chan_id, method='nearest')

Let’s also find the stimulus presentations in this window

stims_in_window = stim_presentations[(stim_presentations.start_time>start_time)&(stim_presentations.start_time<end_time) &
                                    (stim_presentations.omitted==False)]
stim_times_in_window = stims_in_window.start_time.values

Finally, we can plot the spike times and stim times along with the LFP for this interval:

_ = plt.plot(lfp_data.time, lfp_data)
_ = plt.plot(times_in_range, np.ones(times_in_range.shape)*3e-4, '.r')
_ = plt.xlabel('Time (s)')
_ = plt.ylabel('LFP (V)')

_ = plt.plot(stim_times_in_window, np.ones(stim_times_in_window.size)*4e-4, 'vg')

plt.legend(['LFP', 'spikes', 'stim times'])
    

<matplotlib.legend.Legend at 0x7f6bccad95e0>

../_images/6b69aeb41cadf1a0045f2e6ca99b6a95d2a79341a7d6136a1a785e944004374e.png

This plot shows clear troughs in the LFP associated with every stimulus presentation. We also see that the spiking for our selected unit is highly entrained to the stimulus presentations.

Now let’s calculate a spike triggered average of the LFP using a subset of spikes for our unit of interest and the align_lfp function we defined above:

rng = np.random.default_rng(seed=42) #set seed for deterministic results
spikes_to_use = rng.choice(spike_times, min((spike_times.size, 1000)), replace=False)
spike_triggered_lfp = align_lfp(lfp, np.arange(-0.1, 0.1, 1/1250), spikes_to_use)

Let’s plot this spike-triggered LFP for a region of the probe centered on this unit’s peak channel:

fig, ax = plt.subplots()
im = ax.pcolor(spike_triggered_lfp.time_from_presentation_onset.values, lfp_chan_depths, 
               spike_triggered_lfp.mean(dim='presentation_id').data, shading='auto')

ax.plot(-0.01, peak_probe_position, '>w')
ax.text(-0.015, peak_probe_position, 'peak channel', c='w', va='center', ha='right')
ax.set_ylim([peak_probe_position-300, peak_probe_position+300])
ax.set_xlabel('Time from spike (s)')
ax.set_ylabel('Channel depth')

Text(0, 0.5, 'Channel depth')

../_images/bd0b7b9896f9d385ad3d04216383b9bb67023c2a1db258f3ea81ada02cd7ea86.png

Aligning the LFP to spikes from our example unit reveals a low resolution spike waveform (note the thin dark stripe at the center of the plot). As an added exercise, you can compare this plot to the high resolution waveform for this unit stored in the session mean_waveforms dictionary. This plot also gives us some information about the network state that was conducive to spiking for this unit. Note the blue LFP depression in the middle of the plot, indicating that this unit tended to spike during periods of relative network excitability. How might this look for other cortical units? What if we limited our analysis to spikes during a particular stimulus or spontaneous activity?

LFP Analysis

Contents

LFP Analysis#

Tutorial overview#

Learning about the available probes#

Loading LFP data#

Aligning LFP data to a stimulus#

Aligning LFP data to units#

Current Source Density#