edeno/jadhav-2016-data-analysis

Code for analysis on data from Jadhav et al. 2016

Eric L. Denovellis

About me

I am a computational research scientist in Dr. Loren Frank's lab at UCSF. My work focuses on developing scalable, interpretable algorithms and tools to decode, categorize and visualize neural representations. This work extends and applies marked point process switching state space models I developed during my postdoc with Dr. Uri T. Eden. I work closely with experimental collaborators to ensure these algorithms and tools are usable on large scale data.

I also build open source software packages for neural data analysis. These include:

• replay_trajectory_classification a Python package for decoding spatial position represented by neural activity and categorizing the type of trajectory.
• Spyglass is a data analysis framework that facilitates the storage, analysis, visualization, and sharing of neuroscience data to support reproducible research. It is designed to be interoperable with the NWB format (a data standard for neurophysiology) and integrates open-source tools such as SpikeInterface and DeepLabCut into a coherent framework.
• spectral_connectivity is a Python software package that computes multitaper spectral estimates and frequency-domain brain connectivity measures such as coherence, spectral granger causality, and the phase lag index using the multitaper Fourier transform.

Previously, I completed my PhD in computational neuroscience at Boston University with Drs. Daniel H. Bullock and Earl K. Miller. There I developed computational tools and models to understand how prefrontal cortex supports the underlying neural comptuations necessary to switch between contexts. Specifically:

• I showed how synchronous network oscillations in the prefrontal cortex provide a mechanism to flexibly coordinate context representations between groups of neurons during task switching
• I used generalized additive models to show that anterior cingulate neurons can represent context, but do not play a significant role in switching between contexts
• Finally, I developed a set of web-enabled interactive visualization tools designed to provide a multi-dimensional integrated view of electrophysiological datasets.

Contact Me: Email

jadhav-2016-data-analysis's Issues

Are the candripples files the consensus ripples from Loren Frank's lab?

Decide whether it is easier to use the Frank Lab Data Filter Framework or if it's easier to code in own format

I would like to make the switch to Python, but the filter framework is written in Matlab. The framework might be hard to parallelize over in the future and it has a bunch of 'eval' statements that make it hard for Matlab to optimize the code. On the other hand, it's already written and I know it'll give me data in the same format as the original paper.

Setup meeting with Loren Frank

Which frequencies should we focus on?

There are a number of frequencies that we could target with our spectral analysis:

• Ripple frequency (150-250 Hz)
• High Gamma (50-100 Hz)
• Low Gamma (30-50 Hz)
• Beta (12-30 Hz)
• Theta/Alpha (4-12 Hz)

Loren thinks that Low Gamma and Theta are of more interest in a coherence analysis. He's worried that increased ripple frequency power simply reflects an increase in spiking from increased input to both areas.

He also suggested that we look at spectrograms in PFC and maybe use that to determine what frequencies to look at.

Compute ripple evoked response in PFC

Check if there's any evoked potential from the ripple event inducing coherence in PFC

Compute point process spectra of the start times of the ripple events

Rhythms of the ripple start times
Are there enough ripples?

Unsure on how to find information about a single trial

Such as whether it was rewarded or not rewarded, inbound or outbound

Compute Loeve spectrum to examine cross-frequency coupling

Conjugate of complex spectrum in one frequency times the complex spectrum of another . Identifies the ability of phase / magnitude in one frequency to predict the other.

Also look at residuals

Unsure how to find the trials within the epochs

Have not seen a data structure that allows you to get the start and stop times for each trial

Scaling up analysis to process multiple days

Eventually we will have to scale up to process all the day. One possibility is using dask. Examples:

• Example1
• Example2

We should use their computed ripples in the spectral analysis for comparison

Email to figure out how novel is each track exposure

Examine phase locking pfc neurons to hippocampal theta

Compute canonical coherence

Better model for reverse time trajectory for the ripple decoding

Xinyi's model assumes that outbound reverse movements are the same as inbound forward movements in the state transition matrix. We could improve this by computing the reverse time trajectory.

Compute Gamma, Theta relative to ripples

Better function documentation

• data_processing.py
• ripple_detection.py
• spectral.py
• ripple_decoding.py
• analysis.py
• test_data_processing.py
• test_ripple_detection.py
• test_spectral.py
• test_ripple_decoding.py
• test_analysis.py

Should we include iCA1 when defining a ripple?

Compute spectrogram of ripples in ripple band for CA1 tetrodes

Sanity check on detection method

Email paper about group delay

Write group delay test code

Rename data_filter module to data_processing

I think data_processing is probably a better name for a module that involves extracting data from their data structures and reshaping it as necessary

Should we use the Frank lab ripples or the Multitaper ripples

Compute Group Delay

Compute coherence phase vs frequency and look a linear component.

Having trouble finding all the LFP recordings

In the EEG folder, there are several types of files:

• files that correspond to the raw LFP recordings ({animal}eeg{day}-{epoch}-{tetrode}).mat
• files that correspond to the ground LFP recordings ({animal}eeggrnd{day}-{epoch}-{tetrode}.mat
• files that correspond to the processed LFP in different frequency bands (e.g. {animal}delta{day}-{epoch}-{tetrode}

But there doesn't seem to be a file for each tetrode. I'm not sure why.

To further complicate things, there is an OriEEG folder, which I assumed meant "Original EEG" files, but seems to only include tetrodes 17 and above. These seem to be only PFC electrodes (at least for the first animal).

I'm going to try to:

1. See if the Filter framework code gives any clues
2. Maybe email Demetris from Loren's lab

Quantify behavior

The animals start with a novel track and then learn over days. Need to characterize:

• time between reward wells for each session
• amount rewarded for each session
• how these change relative to inbound / outbound trials and over sessions/days

Should expect time between reward wells to get faster over sessions/days and amount rewarded to increase over sessions/days. Inbound trials should be faster than outbound trials.

Use clusterless decoding or sorted spikes?

Clusterless decoding might give us more spikes per ripple to work with (because it includes spikes that couldn't be assigned to a particular neuron)

Criterion for neuron exclusion

Currently, I am excluding all CA1 neurons with zero spikes. The problem is that for low spiking neurons, the training data (times when the rat is running > 4 cm/s) itself might not have spikes.

Compute peri-event time histogram of neurons relative to ripple event

Prepare slides for meeting with Loren Frank

1. High frequency difference in coherence
2. Low frequency difference in coherence
3. Ripple triggered averages
4. Group delay?

Is the encoding model fitting properly?

Use occupancy normalized histograms overlayed on model fit
Maybe simulate data from the model and compare to the real data
Maybe compare to Xinyi’s data

Add session time categorical to the ripple_info data structure

So we can compare ripples at the beginning of the session to ripples at the end of the session. Related to issue #47.

Subtract mean when computing the spectrum when computing power differences / coherence differences

Normalize the ripples before computing the ripple triggered average LFP

See if that affects the ripple triggered average LFP

Write test code for ripple segmenting and spectral functions

Reproduce analysis from Jadhav paper

In order to make sure I understand the data formatting and analysis in the original paper, need to reproduce at least one analysis from the Jadhav paper. Two candidate analyses are:

• Figure 3 - show PFC activity locked to the time of the short wave ripple
• Figure 6 - the GLM decoding of PFC spiking using the CA1 ensemble

Add 1% of firing rate as an offset as a baseline for the conditional intensity

Keeps a spike from being so informative that it forces the decoding filter to immediately choose a state and stay there

Categorize ripple types

• Inbound trials vs. Outbound trials
• Prospective vs. Retrospective Replay
• Right vs Left turn
• Local vs. Remote environment (Remote = some other environment the animal has previously been in)

Take spectrogram of PFC tetrodes on entire experiment

Look for spectral features of PFC LFPs

Add to ripple summary statistics: number of spikes per ripple, number of unique spikes per ripple

Gives us a better idea of the information content of each ripple and our ability to decode

Make an interactive visualization of the ripple decoding

Could be useful for checking the quality of decodes for a large dataset and for demonstration of the algorithm itself.

Figure out how they do their referencing

There are tetrodes labeled CA1Ref, PFCRef, but it's not clear if they are using both references, one reference, or if this is just a potential marker for reference electrode.

Use ripple detection method of Kay 2016?

Our version of their code used their description of ripple detection in
Karlsson, M.P., Frank, L.M., 2009. Awake replay of remote experiences in the hippocampus. Nature Neuroscience 12, 913–918. doi:10.1038/nn.2344:

SWRs were identified on the basis of peaks in the LFP recorded from one channel from each tetrode in the CA3 and CA1 cell layers. The raw LFP data were bandpass-filtered between 150–250 Hz, and the SWR envelope was determined using a Hilbert transform. The envelope was smoothed with a Gaussian (4-ms s.d.). We initially identified SWR events as sets of times when the smoothed envelope stayed above 3 s.d. of the mean for at least 15 ms on at least one tetrode. We defined the entire SWR as including times immediately before and after that threshold crossing event during which the envelope exceeded the mean. Overlapping SWRs were combined across tetrodes, so many events extended beyond the SWR seen on a single tetrode.

Which is what Jadhav 2016 uses as well. But a different detection method is presented in Kay, K., Sosa, M., Chung, J.E., Karlsson, M.P., Larkin, M.C., Frank, L.M., 2016. A hippocampal network for spatial coding during immobility and sleep. Nature 531, 185–190. doi:10.1038/nature17144:

Detection of SWRs was prerequisite for all data analysed in this study, and was performed only when at least three CA1 cell layer recordings were available. Offline, a multisite average approach was used to detect SWRs58. Specifically, LFPs from all available CA1 cell layer tetrodes were filtered between 150–250 Hz, then squared and summed across tetrodes. This sum was smoothed with a Gaussian kernel (σ = 4 ms) and the square root of the smoothed sum was analysed. SWRs were detected when the signal exceeded 2 s.d. of the recording epoch mean for at least 15 ms. SWR periods were then defined as the periods, containing the times of threshold crossing, in which the power trace exceeded the mean. SWR onset was defined as the start of a SWR period. Detection of SWRs was performed only when subjects’ head speed was <4 cm s−1. For SWR-triggered spike raster plots and PSTH plots, a 0.5 s exclusion period was imposed to isolate SWRs occurring only after non-SWR periods; otherwise, analyses of SWRs included all detected SWRs.

DIO pulses are corrupted for a couple of animals

The DIO file has six channels corresponding to motion detection (the animal is near the reward well) or reward delivery for each of the three reward wells located at the end of each arm. So the animals' position (as marked by infrared diodes on top of the microdrive preamp) when these channels are active should correspond to the end of the arms where the reward well is located. However, this is not always the case:

Blue lines indicate the rat's position in the W-track. Red dots correspond to the rat's position when either motion detection or reward delivery is triggered. Data is from rat 'HPa', day 8, epoch 2. See python notebook UnderstandingDIO for more plots and code.

I emailed Demetris and he confirmed there might be an issue:

Unfortunately you will have to ask shantanu about the other issues as I seem to remember him saying something about his dio pules being corrupted for a couple of the animals leading to the displaced position/time issue that you're seeing (so don't worry you're not crazy :)). By now, he may have manually corrected these pulses or reconstructed the dio-like info from the position or something.

I will email Shantanu and see if he has any solutions

Compute coherence by ripple type

• spiking vs no spiking replay events coherence
• first half of the experiment ripples vs second half ripples
• forward vs. reverse
• inbound vs. outbound

Compute spectral granger causality

See:

Geweke, J. 1982 Measurement of linear dependence and feedback between multiple time series. Journal of the American Statistical Association 77, 304-313.

Are there files in the current dataset we can use for marks

Rename segments to peri-event window

More informative name

Apply spectral analysis to different animal / epoch

Add speed criterion to ripple detection