Date: 2016-08-26

The purpose of this notebook is to develop code for exploratory analysis of the LFPs. The goals are to:

  • Plot LFP vs. time
  • Plot the periodograms of the LFPs
  • Plot the spectrograms (power vs. frequency in time) using multitaper spectral analysis
  • Plot the power at a frequency of interest vs. time
  • Plot the standardized power at a frequency of interest vs. time
  • Develop code that detects sharp wave ripples as in Jadhav et al. 2016

In [1]:
%matplotlib inline
%reload_ext autoreload
%autoreload 2
%qtconsole
import os
import sys
import collections
import scipy.io
import scipy.fftpack
import scipy.signal
import scipy.ndimage
import numpy as np
import nitime.algorithms as tsa
import nitime.utils as utils
import nitime.viz
import matplotlib.pyplot as plt
import matplotlib.colors as clr
import seaborn as sns
import pandas as pd
import tqdm

sys.path.append('../src/')
import data_filter as df

Accessing the tetrode data

The following is some code to allow for accessing the needed tetrode data. For this exploratory analysis, I want to focus on the last day for animal HPa when the animal is running on the w-track. First check if we can get the correct epoch info.


In [2]:
Animal = collections.namedtuple('Animal', {'directory', 'short_name'})
num_days = 8
days = range(1, num_days + 1)
animals = {'HPa': Animal(directory='HPa_direct', short_name='HPa')}

epoch_info = df.make_epochs_dataframe(animals, days)
epoch_info.head()


Out[2]:
environment type
animal day epoch_ind
HPa 1 1 presleep sleep
2 lin run
3 NaN rest
4 wtr1 run
5 NaN rest

Get the epoch dataframe corresponding to day 8, w-track


In [3]:
(epoch_info
    .loc[(['HPa'], [8]), :]
    .loc[epoch_info.environment == 'wtr1'])


Out[3]:
environment type
animal day epoch_ind
HPa 8 2 wtr1 run

Next check if we can export the row index information as a tuple so we can use it as a dictionary key to access the tetrode information.


In [4]:
df.get_dataframe_index(epoch_info
    .loc[(['HPa'], [8]), :]
    .loc[epoch_info.environment == 'wtr1'])


Out[4]:
[('HPa', 8, 2)]

Now get the information about the tetrodes. Notice that this structure is a dictionary with keys corresponding to (animal, day, epoch_index) tuples. These keys can access the corresponding dataframes with the tetrode information.


In [5]:
tetrode_info = df.make_tetrode_dataframe(animals)
tetrode_info.keys()


Out[5]:
dict_keys([('HPa', 8, 4), ('HPa', 7, 5), ('HPa', 1, 2), ('HPa', 5, 3), ('HPa', 2, 3), ('HPa', 7, 4), ('HPa', 1, 3), ('HPa', 1, 1), ('HPa', 4, 3), ('HPa', 2, 2), ('HPa', 8, 5), ('HPa', 3, 2), ('HPa', 7, 1), ('HPa', 6, 2), ('HPa', 4, 2), ('HPa', 2, 1), ('HPa', 3, 3), ('HPa', 8, 1), ('HPa', 6, 4), ('HPa', 3, 4), ('HPa', 3, 5), ('HPa', 5, 4), ('HPa', 4, 5), ('HPa', 1, 7), ('HPa', 6, 3), ('HPa', 6, 5), ('HPa', 4, 4), ('HPa', 5, 2), ('HPa', 1, 4), ('HPa', 7, 2), ('HPa', 8, 2), ('HPa', 6, 1), ('HPa', 8, 3), ('HPa', 1, 5), ('HPa', 1, 6), ('HPa', 5, 5), ('HPa', 3, 1), ('HPa', 2, 5), ('HPa', 5, 1), ('HPa', 7, 3), ('HPa', 4, 1), ('HPa', 2, 4)])

Get the tetrode dataframes corresponding to the epoch we want and concatenate them to make a dataframe. In this case, there's only one key so the concatenation step doesn't do anything.


In [6]:
epoch_keys = df.get_dataframe_index(epoch_info
    .loc[(['HPa'], [8]), :]
    .loc[epoch_info.environment == 'wtr1'])

cur_tetrode_info = pd.concat([tetrode_info[key] for key in epoch_keys])
cur_tetrode_info


Out[6]:
area depth descrip numcells
animal day epoch_ind tetrode_number
HPa 8 2 1 CA1 113 riptet 12
2 CA1 121 NaN 0
3 CA1 90 CA1Ref 0
4 CA1 116 riptet 15
5 CA1 116 riptet 0
6 CA1 110 riptet 0
7 CA1 114 riptet 0
8 iCA1 114 riptet 0
9 iCA1 100 riptet 0
10 iCA1 96 NaN 0
11 iCA1 106 riptet 0
12 iCA1 114 riptet 3
13 iCA1 120 NaN 0
14 iCA1 105 riptet 6
15 PFC 93 NaN 0
16 PFC 90 NaN 0
17 PFC 90 NaN 6
18 PFC 90 NaN 0
19 PFC 130 NaN 0
20 PFC 109 NaN 0

Now we have the tetrode information, we can extract keys that correspond to each tetrode. These can be used to access the tetrode LFP data.


In [7]:
tetrode_index = df.get_dataframe_index(cur_tetrode_info)
tetrode_index


Out[7]:
[('HPa', 8, 2, 1),
 ('HPa', 8, 2, 2),
 ('HPa', 8, 2, 3),
 ('HPa', 8, 2, 4),
 ('HPa', 8, 2, 5),
 ('HPa', 8, 2, 6),
 ('HPa', 8, 2, 7),
 ('HPa', 8, 2, 8),
 ('HPa', 8, 2, 9),
 ('HPa', 8, 2, 10),
 ('HPa', 8, 2, 11),
 ('HPa', 8, 2, 12),
 ('HPa', 8, 2, 13),
 ('HPa', 8, 2, 14),
 ('HPa', 8, 2, 15),
 ('HPa', 8, 2, 16),
 ('HPa', 8, 2, 17),
 ('HPa', 8, 2, 18),
 ('HPa', 8, 2, 19),
 ('HPa', 8, 2, 20)]

LFP - time domain

Now retrieve the tetrode data and plot the first tetrode to make sure we have the right data.


In [8]:
def tetrode_title(tetrode_index_tuple, cur_tetrode_info):
    return '{brain_area} Tetrode #{tetrode_number}'.format(tetrode_number=tetrode_index_tuple[-1],
                     brain_area=cur_tetrode_info.loc[tetrode_index_tuple]['area'])

lfp_data = df.get_LFP_data(tetrode_index, animals)
lfp_data[0].plot(title=tetrode_title(tetrode_index[0], cur_tetrode_info), figsize=(12,3))


Out[8]:
<matplotlib.axes._subplots.AxesSubplot at 0x114da96a0>

Now plot all the tetrodes


In [9]:
num_rows = int(np.ceil(len(tetrode_index) / 7))
fig, axis_handles = plt.subplots(num_rows, 7, figsize=(21, 9), sharex=True, sharey=True)
for tetrode_ind, axis_handle in enumerate(axis_handles.flatten()):
    try:
        lfp_data[tetrode_ind].plot(ax=axis_handle,
                                   title=tetrode_title(tetrode_index[tetrode_ind], cur_tetrode_info),
                                  )
    except IndexError:
        axis_handle.xaxis.set_visible(False)
        
plt.suptitle('LFP Time Series: HPa, Day #8, Epoch #2', fontsize=24)


Out[9]:
<matplotlib.text.Text at 0x114c6bba8>

Periodograms

The following is hand coded to make sure I understand the calculation for the periodogram.


In [10]:
def plot_spectrum_by_hand(data_frame, sampling_frequency, axis_handle):
    
    data = data_frame['electric_potential'] - data_frame['electric_potential'].mean()
    number_samples = data.size
    number_fft_bins = np.floor(number_samples / 2) + 1
    nyquist = 0.5 * sampling_frequency

    dft = scipy.fftpack.fft(data)
    frequencies = np.linspace(0.0, nyquist, number_fft_bins)
    power = np.abs(dft[0:number_fft_bins])**2 / (number_fft_bins * sampling_frequency)
    axis_handle.plot(frequencies, power)
    plt.xlim((0, 200))
    plt.xticks([0, 6, 12, 20, 50, 90, 150, 200])

fig, axis_handle = plt.subplots(1, 1, figsize=(12, 9))
sampling_frequency = 1500
plot_spectrum_by_hand(lfp_data[0], sampling_frequency, axis_handle)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power (dB)')


/Users/edeno/anaconda3/envs/Jadhav-2016-Data-Analysis/lib/python3.5/site-packages/ipykernel/__main__.py:10: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
Out[10]:
<matplotlib.text.Text at 0x1187d71d0>

Compare to scipy.signal package periodogram


In [11]:
def plot_spectrum_scipy(data_frame, sampling_frequency, axis_handle):
    frequencies, power_spectral_density = scipy.signal.periodogram(data_frame['electric_potential'], 
                                                                   fs=sampling_frequency,
                                                                   nfft=data_frame['electric_potential'].size,
                                                                   detrend='constant')
    axis_handle.plot(frequencies, power_spectral_density)
    plt.xlim([0, 200])
    plt.xticks([0, 6, 12, 20, 50, 90, 150, 200])
    
fig, axis_handle = plt.subplots(1, 1, figsize=(12, 9))
sampling_frequency = 1500
plot_spectrum_scipy(lfp_data[0], sampling_frequency, axis_handle)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power')


Out[11]:
<matplotlib.text.Text at 0x11765a518>

Compare to the nitime package periodogram


In [12]:
def plot_spectrum_nitime(data_frame, sampling_frequency, axis_handle):
    frequencies, power_spectral_density = tsa.periodogram(data_frame['electric_potential'] - data_frame['electric_potential'].mean(),
                                                          Fs=sampling_frequency,
                                                          normalize=True)
    axis_handle.plot(frequencies, power_spectral_density)
    
fig, axis_handle = plt.subplots(1, 1, figsize=(12, 9))
sampling_frequency = 1500
plot_spectrum_nitime(lfp_data[0], sampling_frequency, axis_handle)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power')
plt.xlim([0, 200])
plt.xticks([0, 6, 12, 20, 50, 90, 150, 200])


Out[12]:
([<matplotlib.axis.XTick at 0x15c7c52b0>,
  <matplotlib.axis.XTick at 0x1187bc828>,
  <matplotlib.axis.XTick at 0x15c7ce860>,
  <matplotlib.axis.XTick at 0x15c7fe320>,
  <matplotlib.axis.XTick at 0x15c7fee10>,
  <matplotlib.axis.XTick at 0x145f51940>,
  <matplotlib.axis.XTick at 0x145f4e470>,
  <matplotlib.axis.XTick at 0x145f4ef60>],
 <a list of 8 Text xticklabel objects>)

Now plot all the periodgrams.


In [13]:
num_rows = int(np.ceil(len(tetrode_index) / 7))
fig, axis_handles = plt.subplots(num_rows, 7, figsize=(21, 9), sharex=True, sharey=True)
for tetrode_ind, axis_handle in enumerate(axis_handles.flatten()):
    try:
        plot_spectrum_nitime(lfp_data[tetrode_ind], sampling_frequency, axis_handle)
        axis_handle.set_title(tetrode_title(tetrode_index[tetrode_ind], cur_tetrode_info))
    except IndexError:
        axis_handle.xaxis.set_visible(False)
        
plt.suptitle('LFP Periodograms: HPa, Day #8, Epoch #2', fontsize=24)


Out[13]:
<matplotlib.text.Text at 0x15c7f3dd8>

Multitaper Spectral Analysis

Let's try to get a better estimate with the multitaper spectrum. Using tapers gives better control over the amount of smoothing in the frequency domain, reducing the bias of the estimate of the signal. The tapers can be designed to tradeoff between the narrow-band bias (the frequency resolution) and the broad-band bias (the spectral leakage from the sidelobes of the smoothing kernel). Using multiple orthogonal tapers that overlap in time-frequency subdomains and averaging over the tapers reduces the variance.

Multitaper spectrum analysis requires three parameters to be set (the frequency resolution, the time window, and the number of tapers). Smaller frequency resolutions are needed to resolve lower frequency rhythms. Smaller time windows are needed for quickly changing non-stationary processes. More tapers are needed to give better control over frequency domain smoothing and reduce the variance of the estimate.

For hippocampal data, the important rhythms are the theta (6-12 Hz), the low gamma (20-50 Hz), the high gamma (50-90 Hz), and the sharp wave ripple band (150-250 Hz).

For prefrontal cortex, the important rhythms are alpha (8-12 Hz), beta (15-30 Hz), and gamma (30-60 Hz).

For the alpha/theta bands we need a good frequency resolution like 2 Hz. This implies a half-bandwidth ($W$) of 1 Hz. The time-half-bandwidth product ($T * W$) has to be greater than 1 ($T*W \ge 1$) where $T$ is the length of time of the signal, so $T = 1$ second. The optimal number of tapers is $2 * T * W - 1$, so we need one taper.

For the beta/gamma bands, we can reduce the frequency resolution to 5 or 10 Hz. For the sharp wave ripple band, we could decrease this to 25 or 50 Hz.

The following table lists the different tradeoffs for frequency resolution, window time length and number of tapers:

Frequency Resolution (Hz) Window Time Length (sec) Number of Tapers
2 1.000 1
5 0.500 1
10 0.200 1
25 0.080 1
50 0.040 1
2 3.000 5
5 1.200 5
10 0.600 5
25 0.240 5
50 0.120 5

It would probably be best to run the multitaper analysis with a couple of different sets of parameters in order to capture the different aspects of the signal. For the alpha/theta, use 2 Hz frequency resolution, 1.000 second window, 1 taper. For beta/gamma/sharp wave ripple, use 10 Hz frequency resolution, 0.600 second window, 5 tapers. For just the timing of the sharp wave ripple, 50 Hz frequency resolution, 0.040 second window, 1 taper.

Prototyping code for the multitaper power spectral density

The following is just some prototype code for computing the multitaper power spectral density over successive windows. The plot corresponds to the power spectral density for each window.

Some notes on the code:

The time domain signal is centered at zero for each window.

The sliding windows share only the last time point (due to how the Pandas package dataframe selection includes both the start and end points given. This may have to change in the future. Alternatively, could put in a paramter that allows for more sliding window overlap.

Lastly, the last time window is excluded because there aren't enough data points. If more window overlap was allowed, this would increase use of the data. Not sure if there's a better solution.


In [14]:
time_window_start = lfp_data[0].index[0]
time_window_duration = 1.00

while time_window_start + time_window_duration < lfp_data[0].index[-1]:
    try:
        windowed_data_frame = lfp_data[0].loc[time_window_start:(time_window_start + time_window_duration), :]
        frequencies, power_spectral_density, _ = tsa.multi_taper_psd(windowed_data_frame['electric_potential'] - windowed_data_frame['electric_potential'].mean(),
                                                                      BW=2,
                                                                      Fs=sampling_frequency,
                                                                      adaptive=False,
                                                                      jackknife=False
                                                                      )
        plt.plot(frequencies, power_spectral_density)
        plt.xlim([0, 20])
        plt.xticks([0, 6, 12, 20])
        plt.ylabel('Power Spectral Density (dB)')
        plt.xlabel('Frequency (Hz)')
        time_window_start += time_window_duration
    except ValueError:
        pass


Now take the prototype code, turn it into a generator function that returns a power spectral density dataframe for each time window. The dataframe consists of the index (time, frequency), where time corresponds to the mid-point time of the window. The frequencies are the frequencies returned from the nitime multi-taper psd function. The dataframe has one column corresponding to power, which is the power spectral density calculated the nitime package.

The dataframe generator can be looped through and concatenated to get the spectrogram dataframe.


In [15]:
def make_timefrequency_dataframe(lfp_data_frame, time_window_duration, time_window_step, frequency_resolution, sampling_frequency):
    ''' Generator function that returns a power spectral density data frame for each time window
    '''
    time_window_start = 0
    number_points_time_window = int(np.fix(time_window_duration * sampling_frequency))
    number_points_time_step = int(np.fix(time_window_step * sampling_frequency))

    while time_window_start + number_points_time_window < len(lfp_data_frame):
        try:
            time_window_end = time_window_start + number_points_time_window
            windowed_data_frame = lfp_data_frame.iloc[time_window_start:time_window_end, :]
            frequencies, power_spectral_density, _ = tsa.multi_taper_psd(windowed_data_frame['electric_potential'] - windowed_data_frame['electric_potential'].mean(),
                                                                          BW=frequency_resolution,
                                                                          Fs=sampling_frequency,
                                                                          adaptive=False,
                                                                          jackknife=False
                                                                          )
            yield (pd.DataFrame({
                                'power': power_spectral_density,
                                'frequency': np.round(frequencies, decimals=3)
                                })
                                .assign(time=lambda x: np.round(lfp_data_frame.index.values[time_window_start] + (time_window_duration / 2), decimals=4))
                                )
            time_window_start += number_points_time_step
        except ValueError:
            # Not enough data points
            pass

def get_spectrogram_dataframe(lfp_data_frame, sampling_frequency=1000, frequency_resolution=2, time_window_duration=1, time_window_step=0):
    ''' Returns a pandas dataframe with the information for a spectrogram. Sampling frequency and frequency
    resolution inputs are given in Hertz. Time window duration and steps are given in seconds.
    '''
    if time_window_step == 0:
        time_window_step = time_window_duration
    return pd.concat((time_window
                      for time_window in make_timefrequency_dataframe(lfp_data_frame,
                                                                      time_window_duration,
                                                                      time_window_step,
                                                                      frequency_resolution,
                                                                      sampling_frequency)))

In [16]:
frequency_resolution = 2
time_window_duration = 1
sampling_frequency = 1500
time_window_step = time_window_duration / 2
spectrogram = get_spectrogram_dataframe(lfp_data[0],
                                        frequency_resolution=frequency_resolution,
                                        time_window_duration=time_window_duration,
                                        sampling_frequency=sampling_frequency,
                                        time_window_step=time_window_step)

spectrogram.head()


Out[16]:
frequency power time
0 0.0 9.745951 2713.4948
1 1.0 111.634740 2713.4948
2 2.0 346.019641 2713.4948
3 3.0 562.672609 2713.4948
4 4.0 19.718734 2713.4948

If we pivot the spectrogram dataframe, we get a two-dimensional array with frequency on the y-axis and time on the x-axis


In [17]:
spectrogram.pivot('frequency', 'time', 'power')


Out[17]:
time 2713.4948 2713.9948 2714.4948 2714.9948 2715.4948 2715.9948 2716.4948 2716.9948 2717.4948 2717.9948 ... 3914.9949 3915.4949 3915.9949 3916.4949 3916.9949 3917.4949 3917.9949 3918.4949 3918.9949 3919.4949
frequency
0.0 9.745951 6.362539 27.732212 5.917123 3.650063 6.940113 2.061648e+00 120.440799 0.075193 112.697688 ... 64.974900 69.496085 13.202142 59.322530 3.377219 1.688537 1.329371 71.483425 97.821297 16.453320
1.0 111.634740 33.151833 108.732378 46.514179 284.128803 391.463048 5.565351e+02 1358.339987 479.811706 1377.505036 ... 948.693800 925.922444 503.737218 438.462336 20.390079 118.975513 86.669778 503.010180 666.548472 106.570633
2.0 346.019641 848.785288 12.738344 9.481091 688.665054 851.619984 2.959881e+02 350.492282 1304.749580 330.814450 ... 988.348157 508.349305 836.706119 134.971445 268.301811 538.099463 295.807257 614.028046 365.651909 186.510312
3.0 562.672609 314.992579 466.181980 254.743371 619.143095 102.239638 3.094345e+02 228.302585 29.906341 630.482113 ... 495.039078 128.772579 3.674497 89.051187 260.434143 391.007002 174.420333 197.420525 392.761581 846.537442
4.0 19.718734 86.598309 543.763036 98.040628 6.103205 46.264467 1.203846e+02 276.879102 77.677463 141.674710 ... 77.402016 22.178221 76.176241 146.492076 89.846631 41.833086 150.700662 110.447949 423.168836 0.011780
5.0 397.216604 102.756020 1400.722811 1152.485271 462.074643 199.609446 1.081479e+02 48.597044 61.567197 299.046498 ... 41.471532 66.848224 189.433861 380.967584 2.913539 128.175006 73.985897 220.872535 345.723742 352.477938
6.0 1329.227906 652.838610 1676.203324 1056.922428 710.445028 94.315713 1.674088e+02 249.884752 53.954529 150.905760 ... 572.144177 113.987375 45.547592 19.840563 46.080856 246.160359 153.077768 862.139330 1213.893361 155.526944
7.0 149.775340 324.089474 259.459660 50.964451 1.960415 107.644031 1.298348e+01 129.358399 241.066790 649.347790 ... 477.155587 188.004681 86.889048 310.949836 46.986761 40.496541 216.995712 393.011196 79.580547 56.665575
8.0 180.283287 676.374234 62.182004 74.560905 64.306257 2.283609 1.803855e+01 37.733601 18.776231 401.706054 ... 193.482565 527.984673 1366.115527 290.571764 204.744926 104.401611 90.685429 268.352777 11.407603 18.282818
9.0 19.795969 27.483157 80.891232 124.377429 142.461649 65.043676 1.261262e+02 7.318994 363.863635 73.113936 ... 85.823460 213.565066 426.256616 450.823431 156.880429 9.941648 58.327670 32.796785 13.094284 219.244265
10.0 52.566635 30.542490 35.560144 50.907878 27.069773 73.042793 1.923286e+01 77.460200 372.419042 24.435487 ... 75.909619 61.364660 26.580014 116.810209 68.588863 296.527749 2.023204 136.359480 21.363701 272.510683
11.0 5.557012 375.970407 67.708515 96.402872 1.587404 50.384393 2.059075e+01 0.454116 315.433911 70.316006 ... 21.967053 41.975166 138.634615 160.356296 119.062571 207.829779 86.711521 336.916148 99.592873 175.369534
12.0 99.237672 527.950253 452.528407 156.034718 81.811866 17.823252 1.244811e+01 94.585793 100.470609 63.229253 ... 27.039980 49.351843 3.659085 250.392673 122.700904 167.587651 281.040144 39.901230 28.060599 98.378074
13.0 14.769713 12.447772 12.609316 51.500147 32.816278 348.422750 5.227877e+01 73.548709 391.584056 54.856002 ... 119.645033 104.228519 1.699787 102.339377 236.758359 55.366633 21.911145 27.612175 7.272319 261.900553
14.0 98.915065 30.399960 38.176258 81.006823 12.154309 154.753940 6.005461e+02 13.540234 427.998163 717.554805 ... 387.697866 268.618517 257.117436 298.780391 32.765114 2.625921 112.594975 63.053568 48.825564 7.372076
15.0 112.640500 221.391803 149.716117 234.388908 87.509922 74.536194 1.351273e+03 38.964744 142.220535 48.271494 ... 25.223420 27.031929 186.163941 61.324700 232.238821 99.368757 43.727885 125.891168 99.847659 69.163483
16.0 60.719843 12.094028 89.699629 62.512996 71.043388 547.886616 8.066508e+02 137.694936 46.987822 103.618452 ... 128.209268 403.101920 739.368176 576.936031 16.380194 61.887375 106.616435 30.029558 60.550539 11.318492
17.0 239.952717 55.902046 58.517427 12.220604 34.000839 97.595441 3.584163e+00 65.871350 100.024251 39.448441 ... 197.683030 121.439637 82.385711 79.301183 117.163878 111.014618 46.358928 82.910269 91.081339 34.954139
18.0 68.167446 105.387951 15.001185 25.320923 26.466702 90.576718 2.400629e+02 234.140682 81.471712 29.421709 ... 14.327301 150.853583 119.646677 16.422462 168.936223 51.615366 234.911157 189.747350 0.924540 57.782928
19.0 14.321445 15.361808 32.007144 12.274784 61.465770 57.135253 1.809940e+02 7.932123 44.364090 116.544359 ... 178.845470 119.267900 36.465803 4.412988 54.580077 38.151891 54.884777 71.294483 190.656369 31.929063
20.0 19.176213 48.486541 184.747947 17.793353 16.679643 51.952033 3.247871e+01 55.510795 22.940566 24.400632 ... 210.883563 62.868535 11.527567 12.659713 14.116017 64.446798 57.156752 63.344325 71.762560 1.461584
21.0 154.169214 10.359935 82.898778 119.611765 150.437606 13.905484 7.885776e-02 56.531522 20.043479 24.037399 ... 173.473596 2.177004 23.414304 70.674929 8.337117 128.459175 20.283298 15.536589 33.679176 27.014513
22.0 110.609571 0.508947 72.955577 413.302498 53.357926 3.838297 1.838876e+01 37.116878 96.969825 54.915808 ... 47.437781 50.208663 46.481609 117.216907 16.585577 13.986873 55.092485 49.007034 42.948073 92.861317
23.0 28.919507 160.878843 133.231658 0.566630 202.431566 84.350289 7.133918e+01 60.534977 81.117338 1.725753 ... 67.858111 85.592653 24.494443 4.885834 121.121813 43.825715 70.266988 39.838152 164.497875 208.903275
24.0 109.896299 48.817970 539.144586 407.343832 57.999940 85.746850 3.771194e+01 47.602908 41.892671 15.035132 ... 63.961213 72.968853 14.140173 18.884624 26.752403 42.117684 67.102091 72.221112 103.703765 376.378512
25.0 178.918078 35.525015 239.273211 83.308191 11.893576 130.697358 1.602170e+02 125.435436 16.294238 55.569464 ... 9.214343 23.759049 56.909478 27.436823 20.495451 31.681641 80.727281 36.188241 25.121261 94.221841
26.0 28.400421 58.836455 116.083388 20.768503 10.348768 72.844583 5.013077e+01 138.760949 100.288608 8.113784 ... 44.858971 109.740623 34.017051 48.328957 18.823240 36.106081 209.639651 123.205461 33.378453 199.336429
27.0 23.807970 63.946803 18.905584 11.008105 95.977803 170.662570 1.287710e+02 76.143650 77.178339 31.790652 ... 45.120369 26.347127 18.435661 11.250656 105.497715 56.752714 58.044934 4.566944 88.076046 107.237371
28.0 112.101493 85.085504 285.177630 290.548302 5.379892 97.277672 3.542986e+01 77.064569 48.002070 4.589714 ... 23.916514 17.511958 10.556107 5.993122 103.978756 20.665985 4.873379 33.581014 33.127580 119.915234
29.0 331.276779 100.714179 109.224560 126.362506 72.417803 65.565488 1.799818e+01 43.656241 53.802096 10.074241 ... 91.638200 16.172981 31.375761 121.212588 51.674416 4.663049 6.398389 36.181655 122.228745 30.185455
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
721.0 0.013387 0.017908 0.019062 0.018429 0.021883 0.014802 4.713829e-03 0.003496 0.057824 0.026445 ... 0.020356 0.006520 0.006726 0.036257 0.011737 0.007562 0.039269 0.001951 0.005548 0.021972
722.0 0.003841 0.019754 0.005968 0.007936 0.048876 0.013745 7.647453e-03 0.002040 0.042659 0.019363 ... 0.001310 0.000176 0.003855 0.045016 0.019828 0.012914 0.045036 0.005410 0.037495 0.069620
723.0 0.016842 0.002967 0.001716 0.017327 0.017496 0.006654 2.896689e-02 0.008905 0.025854 0.006308 ... 0.000455 0.005510 0.021922 0.003743 0.021725 0.052170 0.006458 0.008332 0.010580 0.005279
724.0 0.000296 0.004692 0.005225 0.014377 0.008603 0.001491 1.544831e-03 0.007468 0.011800 0.004243 ... 0.006675 0.028668 0.025879 0.009804 0.017651 0.060326 0.066576 0.052905 0.013569 0.036082
725.0 0.025287 0.036185 0.004729 0.000507 0.019840 0.007324 2.497724e-03 0.020424 0.020001 0.015869 ... 0.005347 0.020606 0.002544 0.002199 0.035520 0.000471 0.005172 0.001852 0.000340 0.013955
726.0 0.004409 0.040796 0.019555 0.023479 0.003548 0.000098 8.294513e-03 0.005938 0.010381 0.003744 ... 0.000203 0.033487 0.069377 0.007448 0.005437 0.064858 0.052243 0.014227 0.011135 0.002734
727.0 0.004448 0.025114 0.005275 0.021209 0.000265 0.005488 2.327396e-02 0.011181 0.026792 0.015090 ... 0.004374 0.021605 0.263672 0.018040 0.036822 0.010406 0.006155 0.026559 0.036849 0.011159
728.0 0.010482 0.017097 0.005752 0.024737 0.025046 0.009734 4.155744e-03 0.002375 0.015770 0.001158 ... 0.014562 0.019428 0.197505 0.060334 0.052856 0.045472 0.040070 0.018705 0.051510 0.032492
729.0 0.040336 0.014877 0.016939 0.017933 0.037423 0.022281 6.458350e-03 0.010746 0.041120 0.003955 ... 0.012314 0.007226 0.037186 0.003787 0.021620 0.006029 0.035486 0.013439 0.007525 0.013161
730.0 0.003027 0.013817 0.008601 0.036826 0.000876 0.021535 3.889325e-02 0.016752 0.012830 0.001895 ... 0.008044 0.041284 0.000210 0.000128 0.004397 0.012334 0.002168 0.012925 0.012972 0.004189
731.0 0.020549 0.002444 0.004385 0.003116 0.037860 0.001501 1.218815e-02 0.022409 0.038378 0.016148 ... 0.000750 0.007594 0.038015 0.050188 0.006790 0.026791 0.026871 0.013474 0.014026 0.024532
732.0 0.010455 0.000367 0.004854 0.008578 0.026229 0.005818 1.681702e-02 0.001610 0.033311 0.018072 ... 0.016111 0.002352 0.001810 0.028220 0.003377 0.001427 0.000009 0.012773 0.017598 0.048206
733.0 0.002139 0.004034 0.004481 0.001642 0.020969 0.006642 1.106281e-02 0.021119 0.085938 0.015265 ... 0.015276 0.006891 0.002652 0.035124 0.024798 0.033940 0.026460 0.005772 0.007654 0.014940
734.0 0.001931 0.008203 0.003070 0.008024 0.003796 0.000197 5.630554e-03 0.008627 0.027995 0.010951 ... 0.001719 0.023935 0.026655 0.004033 0.008812 0.045392 0.010333 0.047101 0.044815 0.009829
735.0 0.037817 0.031926 0.004296 0.008936 0.027181 0.066413 8.860795e-03 0.023909 0.074975 0.003324 ... 0.006385 0.016568 0.016404 0.003618 0.026813 0.052073 0.003332 0.059642 0.003747 0.006765
736.0 0.030729 0.010253 0.003703 0.006235 0.006792 0.038073 6.916085e-03 0.022086 0.001311 0.029151 ... 0.001450 0.016967 0.021677 0.024114 0.011530 0.001354 0.001064 0.013686 0.004263 0.007608
737.0 0.014218 0.003578 0.012446 0.012195 0.007182 0.002196 5.323248e-03 0.009591 0.007760 0.006379 ... 0.005376 0.021573 0.116578 0.028660 0.030480 0.119683 0.046114 0.000675 0.007735 0.010166
738.0 0.019782 0.011105 0.002132 0.001936 0.019627 0.026371 6.179332e-04 0.009052 0.004434 0.009030 ... 0.002394 0.019649 0.005064 0.010468 0.005751 0.093897 0.112120 0.023591 0.007638 0.002248
739.0 0.022448 0.000718 0.001752 0.010235 0.002081 0.017649 1.436990e-02 0.017277 0.000481 0.006167 ... 0.003145 0.007420 0.026897 0.030176 0.022268 0.015632 0.012235 0.018069 0.031788 0.013129
740.0 0.023309 0.017980 0.007542 0.014859 0.001225 0.005065 1.692192e-02 0.043630 0.091556 0.012198 ... 0.005341 0.006462 0.006531 0.000560 0.023984 0.018740 0.020809 0.018234 0.035940 0.020403
741.0 0.008243 0.026900 0.001122 0.009487 0.021145 0.042928 1.400085e-02 0.045522 0.013962 0.004500 ... 0.006504 0.011556 0.015679 0.005680 0.024800 0.048459 0.042263 0.017581 0.002283 0.007735
742.0 0.021256 0.027198 0.004546 0.024359 0.003398 0.028504 6.908547e-03 0.001271 0.017061 0.069366 ... 0.000755 0.010460 0.020026 0.015665 0.011719 0.073728 0.008961 0.051900 0.015426 0.005036
743.0 0.015359 0.025228 0.006196 0.013017 0.005327 0.037446 4.433076e-03 0.000387 0.014467 0.047387 ... 0.001313 0.001511 0.127608 0.038969 0.026640 0.006300 0.040943 0.014178 0.014470 0.013243
744.0 0.036200 0.002806 0.005529 0.009226 0.081862 0.009620 9.933254e-03 0.011251 0.016695 0.006353 ... 0.012125 0.002302 0.135406 0.042447 0.065981 0.002812 0.016367 0.001262 0.010590 0.020436
745.0 0.011485 0.012421 0.054459 0.044448 0.007674 0.009439 2.816797e-03 0.020548 0.086613 0.009439 ... 0.002092 0.027317 0.101372 0.129359 0.014072 0.026015 0.036062 0.011149 0.012151 0.033857
746.0 0.002467 0.011151 0.012886 0.014325 0.004955 0.024054 2.469937e-02 0.039108 0.001170 0.054156 ... 0.003340 0.010389 0.053624 0.004329 0.006207 0.004824 0.017364 0.044851 0.015284 0.015272
747.0 0.006959 0.011338 0.010493 0.015155 0.016333 0.001366 6.499052e-02 0.003990 0.041532 0.006024 ... 0.018896 0.013039 0.001967 0.008680 0.013849 0.024358 0.033802 0.013081 0.030855 0.007428
748.0 0.008008 0.016889 0.001006 0.007685 0.014317 0.005912 2.688338e-02 0.014295 0.017466 0.003488 ... 0.015946 0.003108 0.025067 0.000033 0.045624 0.006203 0.031427 0.015639 0.001520 0.003088
749.0 0.043917 0.030995 0.007167 0.000574 0.000337 0.000304 3.094592e-02 0.051671 0.004806 0.020852 ... 0.010823 0.008318 0.020401 0.027351 0.076724 0.041111 0.050486 0.014816 0.016197 0.012735
750.0 0.002929 0.000544 0.000182 0.002678 0.000058 0.001083 3.562833e-08 0.001992 0.021237 0.000784 ... 0.008277 0.031355 0.028355 0.002859 0.000469 0.073059 0.001487 0.002195 0.005977 0.018130

751 rows × 2413 columns


In [18]:
def _get_time_freq_from_spectrogram(spectrogram_dataframe):
    return (np.unique(spectrogram_dataframe['time']), np.unique(spectrogram_dataframe['frequency']))

def plot_spectrogram(spectrogram_dataframe, axis_handle):
    time, freq = _get_time_freq_from_spectrogram(spectrogram_dataframe)
    mesh = axis_handle.pcolormesh(time, freq, spectrogram_dataframe.pivot('frequency', 'time', 'power'),
                                  cmap='viridis',
                                  shading='gouraud')
    axis_handle.set_ylabel('Frequency (Hz)')
    axis_handle.set_xlabel('Time (seconds)')
    axis_handle.set_xlim([time.min(), time.max()])
    axis_handle.set_ylim([freq.min(), freq.max()])
    return mesh
    
    
fig, axis_handles = plt.subplots(2, 1, figsize=(15, 12), gridspec_kw={'height_ratios':[1, 3]})

lfp_data[0].plot(title=tetrode_title(tetrode_index[0], cur_tetrode_info), ax=axis_handles[0])

mesh = plot_spectrogram(spectrogram.loc[(spectrogram.frequency <= 20)], axis_handles[1])
yticks = [6, 12, 20]
axis_handles[1].set_yticks(yticks)
fig.colorbar(mesh, label='Power Spectral Density', orientation='horizontal')
plt.tight_layout()



In [19]:
def plot_spectrum(spectrogram_dataframe, axis_handle):
    (spectrogram_dataframe.groupby('frequency')
                         .agg(np.mean)
                         .drop('time', 1)
                         .plot(ax=axis_handle))   

fig, axis_handle = plt.subplots(1, 1, figsize=(12, 9))
plot_spectrum(spectrogram.loc[spectrogram.frequency <= 20], axis_handle)
xticks = [6, 12, 20]
axis_handle.set_xticks(xticks)


Out[19]:
[<matplotlib.axis.XTick at 0x15cbd8978>,
 <matplotlib.axis.XTick at 0x138800d68>,
 <matplotlib.axis.XTick at 0x15e1345f8>]

Now compute for all tetrodes


In [20]:
frequency_resolution = 2
time_window_duration = 1
sampling_frequency = 1500
time_window_step = time_window_duration / 2

spectrograms_1 = [get_spectrogram_dataframe(tetrode,
                           frequency_resolution=frequency_resolution,
                           time_window_duration=time_window_duration,
                           sampling_frequency=sampling_frequency,
                           time_window_step=time_window_step) for tetrode in tqdm.tqdm_notebook(lfp_data)]




In [21]:
# Spectrum
num_rows = int(np.ceil(len(tetrode_index) / 7))
fig, axis_handles = plt.subplots(num_rows, 7, figsize=(21, 9), sharex=True, sharey=True)
for tetrode_ind, axis_handle in enumerate(tqdm.tqdm_notebook(axis_handles.flatten())):
    try:
        selection = spectrograms_1[tetrode_ind].frequency <= 20
        mesh = plot_spectrum(spectrograms_1[tetrode_ind].loc[selection], axis_handle)
        xticks = [6, 12, 20]
        axis_handle.set_xticks((2991, 2992))
        axis_handle.set_title(tetrode_title(tetrode_index[tetrode_ind], cur_tetrode_info))
    except IndexError:
        axis_handle.xaxis.set_visible(False)
        
plt.suptitle('LFP Spectrograms: HPa, Day #8, Epoch #2', fontsize=24)


# Spectrograms
num_rows = int(np.ceil(len(tetrode_index) / 7))
fig, axis_handles = plt.subplots(num_rows, 7, figsize=(21, 9), sharex=True, sharey=True)
for tetrode_ind, axis_handle in enumerate(tqdm.tqdm_notebook(axis_handles.flatten())):
    try:
        selection = spectrograms_1[tetrode_ind].frequency <= 20
        mesh = plot_spectrogram(spectrograms_1[tetrode_ind].loc[selection], axis_handle)
        axis_handle.set_title(tetrode_title(tetrode_index[tetrode_ind], cur_tetrode_info))
        axis_handle.set_xticks((2991, 2992))
    except IndexError:
        axis_handle.xaxis.set_visible(False)
        
plt.suptitle('LFP Spectrum: HPa, Day #8, Epoch #2', fontsize=24)



Out[21]:
<matplotlib.text.Text at 0x15b93db38>

It's hard to tell if the colorscale is wonky on the above plots. Let's try something with a lower frequency resolution to compare.


In [22]:
frequency_resolution = 10
time_window_duration = 0.600
sampling_frequency = 1500
time_window_step = time_window_duration / 2

spectrograms_2 = [get_spectrogram_dataframe(tetrode,
                           frequency_resolution=frequency_resolution,
                           time_window_duration=time_window_duration,
                           sampling_frequency=sampling_frequency,
                           time_window_step=time_window_step) for tetrode in tqdm.tqdm_notebook(lfp_data)]



Comparing the spectra of the different multi-tapers. Blue is lower frequency resolution but more tapers (10 Hz, 5 tapers). Green is higher frequency resolution, but less tapers (2 Hz, 1 taper).


In [23]:
(spectrograms_2[0].groupby('frequency')
                     .agg(np.mean)
                     .drop('time', 1))


Out[23]:
power
frequency
0.000 176.189006
1.667 343.663910
3.333 332.272248
5.000 347.056242
6.667 267.966693
8.333 218.165062
10.000 191.973831
11.667 142.246020
13.333 101.140350
15.000 96.175841
16.667 94.570901
18.333 90.705775
20.000 79.495114
21.667 73.630956
23.333 69.876630
25.000 65.568301
26.667 60.506151
28.333 55.535413
30.000 51.097367
31.667 47.483094
33.333 43.902816
35.000 40.770930
36.667 37.557853
38.333 35.027582
40.000 32.746456
41.667 30.390330
43.333 28.353626
45.000 26.735022
46.667 25.350251
48.333 24.007347
... ...
701.667 0.015179
703.333 0.015164
705.000 0.015206
706.667 0.015131
708.333 0.015191
710.000 0.015078
711.667 0.015001
713.333 0.014940
715.000 0.014808
716.667 0.014875
718.333 0.014928
720.000 0.014867
721.667 0.014933
723.333 0.014797
725.000 0.014845
726.667 0.014834
728.333 0.014734
730.000 0.014755
731.667 0.014898
733.333 0.014807
735.000 0.014705
736.667 0.014779
738.333 0.014704
740.000 0.014617
741.667 0.014588
743.333 0.014639
745.000 0.014687
746.667 0.014619
748.333 0.014589
750.000 0.007307

451 rows × 1 columns


In [24]:
fig, axis_handle = plt.subplots(1, 1, figsize=(12, 9))
plot_spectrum(spectrograms_2[0].loc[spectrograms_2[0].frequency <= 100], axis_handle)
plot_spectrum(spectrograms_1[0].loc[spectrograms_1[0].frequency <= 100], axis_handle)



In [25]:
frequency_resolution = 10
time_window_duration = 0.600
sampling_frequency = 1500
time_window_step = time_window_duration / 2

spectrogram_2 = get_spectrogram_dataframe(lfp_data[0],
                           frequency_resolution=frequency_resolution,
                           time_window_duration=time_window_duration,
                           sampling_frequency=sampling_frequency,
                           time_window_step=time_window_step)

In [26]:
spectrogram_2.head()


Out[26]:
frequency power time
0 0.000 37.325004 2713.2948
1 1.667 161.332339 2713.2948
2 3.333 281.598864 2713.2948
3 5.000 293.425531 2713.2948
4 6.667 287.148589 2713.2948

Compute all of the spectra / spectrograms for the 10 Hz frequency resolution.


In [27]:
# Spectrum
num_rows = int(np.ceil(len(tetrode_index) / 7))
fig, axis_handles = plt.subplots(num_rows, 7, figsize=(21, 9), sharex=True, sharey=True)
for tetrode_ind, axis_handle in enumerate(tqdm.tqdm_notebook(axis_handles.flatten())):
    try:
        mesh = plot_spectrum(spectrograms_2[tetrode_ind].loc[(spectrograms_2[tetrode_ind].frequency >= 20) & (spectrograms_2[tetrode_ind].frequency <= 100)], axis_handle)
        xticks = [20, 30, 50, 60, 90, 150, 200]
        axis_handle.set_xticks(xticks)
        axis_handle.set_title(tetrode_title(tetrode_index[tetrode_ind], cur_tetrode_info))
    except IndexError:
        axis_handle.xaxis.set_visible(False)
        
plt.suptitle('LFP Spectrograms: HPa, Day #8, Epoch #2', fontsize=24)


# Spectrograms
num_rows = int(np.ceil(len(tetrode_index) / 7))
fig, axis_handles = plt.subplots(num_rows, 7, figsize=(21, 9), sharex=True, sharey=True)
for tetrode_ind, axis_handle in enumerate(tqdm.tqdm_notebook(axis_handles.flatten())):
    try:
        mesh = plot_spectrogram(spectrograms_2[tetrode_ind].loc[(spectrograms_2[tetrode_ind].frequency >= 20) & (spectrograms_2[tetrode_ind].frequency <= 100)], axis_handle)
        axis_handle.set_title(tetrode_title(tetrode_index[tetrode_ind], cur_tetrode_info))
    except IndexError:
        axis_handle.xaxis.set_visible(False)
        
plt.suptitle('LFP Spectrum: HPa, Day #8, Epoch #2', fontsize=24)



Out[27]:
<matplotlib.text.Text at 0x15b88fa58>

The spectra/spectrograms above don't quite look right. Maybe the number of tapers is throwing things off? The next cell compares two spectra with frequency resolution of 2 Hz, but different time window lengths of 1 second and 3 seconds (hence different numbers of tapers, 1 vs. 5 respectively). There should be more smoothing in the frequency domain with more tapers.


In [28]:
frequency_resolution = 2
time_window_duration = 3
sampling_frequency = 1500
time_window_step = time_window_duration / 2

spectrograms_3 = [get_spectrogram_dataframe(lfp_data[0],
                           frequency_resolution=frequency_resolution,
                           time_window_duration=time_window_duration,
                           sampling_frequency=sampling_frequency,
                           time_window_step=time_window_step)]

In [29]:
fig, axis_handle = plt.subplots(1, 1, figsize=(12, 9))
plot_spectrum(spectrograms_1[0].loc[spectrograms_1[0].frequency < 200], axis_handle)
plot_spectrum(spectrograms_3[0].loc[spectrograms_3[0].frequency < 200], axis_handle)


Sharp Wave Ripples

For the sharp wave ripples, want low frequency resolution, good time resolution.


In [30]:
frequency_resolution = 100
time_window_duration = 0.020
sampling_frequency = 1500
time_window_step = time_window_duration / 2

spectrograms_swr = [get_spectrogram_dataframe(tetrode,
                           frequency_resolution=frequency_resolution,
                           time_window_duration=time_window_duration,
                           sampling_frequency=sampling_frequency,
                           time_window_step=time_window_step) for tetrode in tqdm.tqdm_notebook(lfp_data)]




In [31]:
fig, axis_handle = plt.subplots(1, 1, figsize=(12, 9))
mesh = plot_spectrogram(spectrograms_swr[0].loc[(spectrograms_swr[0].frequency >= 100) & (spectrograms_swr[0].frequency <= 400)], axis_handle)
plt.colorbar(mesh, label='Power Spectral Density (dB)')

fig, axis_handle = plt.subplots(1, 1, figsize=(12, 9))
plot_spectrum(spectrograms_swr[0].loc[(spectrograms_swr[0].frequency >= 100) & (spectrograms_swr[0].frequency <= 400)], axis_handle)



In [32]:
fig, axis_handles = plt.subplots(2, 1, figsize=(15, 9), sharex=True)

sharp_wave_power = spectrograms_swr[0].loc[(spectrograms_swr[0].frequency == 200)]
time = sharp_wave_power['time']
axis_handles[0].plot(time, sharp_wave_power['power'])
axis_handles[0].set_title('Sharp Wave Power at 200 Hz')

# zscore
def zscore(x):
    return (x - x.mean()) / x.std()

z_sharp_wave_power = zscore(sharp_wave_power['power'])

axis_handles[1].plot(time, z_sharp_wave_power, label='z-scored')
axis_handles[1].set_title('Z-scored Power')
axis_handles[1].set_xlim((time.min(), time.max()))

sharp_wave_power.head()

# Average across electrodes to get a consensus trace

# Threshold at zscore of 3

# Needs to stay above threshold for at least 15 ms on at least one tetrode

# Start and end of ripple are when the zscore exceeds the mean around the event


Out[32]:
frequency power time
4 200.0 2.211998 2713.0048
4 200.0 0.599675 2713.0148
4 200.0 0.648345 2713.0248
4 200.0 0.084342 2713.0348
4 200.0 0.176088 2713.0448

In [33]:
# Spectrum
num_rows = int(np.ceil(len(tetrode_index) / 7))
fig, axis_handles = plt.subplots(num_rows, 7, figsize=(21, 9), sharex=True, sharey=True)
for tetrode_ind, axis_handle in enumerate(tqdm.tqdm_notebook(axis_handles.flatten())):
    try:
        mesh = plot_spectrum(spectrograms_swr[tetrode_ind].loc[(spectrograms_swr[tetrode_ind].frequency >= 100) & (spectrograms_swr[tetrode_ind].frequency <= 300)], axis_handle)
        axis_handle.set_xticks(xticks)
        axis_handle.set_title(tetrode_title(tetrode_index[tetrode_ind], cur_tetrode_info))
    except IndexError:
        axis_handle.xaxis.set_visible(False)
        
plt.suptitle('LFP Spectrograms: HPa, Day #8, Epoch #2', fontsize=24)


# Spectrograms
num_rows = int(np.ceil(len(tetrode_index) / 7))
fig, axis_handles = plt.subplots(num_rows, 7, figsize=(21, 9), sharex=True, sharey=True)
for tetrode_ind, axis_handle in enumerate(tqdm.tqdm_notebook(axis_handles.flatten())):
    try:

        mesh = plot_spectrogram(spectrograms_swr[tetrode_ind].loc[(spectrograms_swr[tetrode_ind].frequency >= 100) & (spectrograms_swr[tetrode_ind].frequency <= 300)], axis_handle)
        axis_handle.set_title(tetrode_title(tetrode_index[tetrode_ind], cur_tetrode_info))
    except IndexError:
        axis_handle.xaxis.set_visible(False)
        
plt.suptitle('LFP Spectrum: HPa, Day #8, Epoch #2', fontsize=24)



Out[33]:
<matplotlib.text.Text at 0x15a6aebe0>

In [ ]: