Ripple Detection¶

Overview¶

Developer Note: if you may make a PR in the future, be sure to copy this notebook, and use the gitignore prefix temp to avoid future conflicts.

This is one notebook in a multi-part series on Spyglass.

• To set up your Spyglass environment and database, see the Setup notebook
• For additional info on DataJoint syntax, including table definitions and inserts, see the Insert Data notebook

Ripple detection depends on a set of LFPs, the parameters used for detection and the speed of the animal. You will need RippleLFPSelection, RippleParameters, and PositionOutput to be populated accordingly.

Imports¶

In [ ]:

import os
import datajoint as dj
import numpy as np

# change to the upper level folder to detect dj_local_conf.json
if os.path.basename(os.getcwd()) == "notebooks":
    os.chdir("..")
dj.config.load("dj_local_conf.json")  # load config for database connection info

import spyglass.common as sgc
import spyglass.position as sgp
import spyglass.lfp as lfp
import spyglass.lfp.analysis.v1 as lfp_analysis
from spyglass.lfp import LFPOutput
import spyglass.lfp.v1 as sglfp
from spyglass.position import PositionOutput
import spyglass.ripple.v1 as sgrip
import spyglass.ripple.v1 as sgr

# ignore datajoint+jupyter async warnings
import warnings

warnings.simplefilter("ignore", category=DeprecationWarning)
warnings.simplefilter("ignore", category=ResourceWarning)

import os import datajoint as dj import numpy as np # change to the upper level folder to detect dj_local_conf.json if os.path.basename(os.getcwd()) == "notebooks": os.chdir("..") dj.config.load("dj_local_conf.json") # load config for database connection info import spyglass.common as sgc import spyglass.position as sgp import spyglass.lfp as lfp import spyglass.lfp.analysis.v1 as lfp_analysis from spyglass.lfp import LFPOutput import spyglass.lfp.v1 as sglfp from spyglass.position import PositionOutput import spyglass.ripple.v1 as sgrip import spyglass.ripple.v1 as sgr # ignore datajoint+jupyter async warnings import warnings warnings.simplefilter("ignore", category=DeprecationWarning) warnings.simplefilter("ignore", category=ResourceWarning)

[2023-07-28 14:45:50,776][INFO]: Connecting root@localhost:3306
[2023-07-28 14:45:50,804][INFO]: Connected root@localhost:3306

Selecting Electrodes¶

First, we'll pick the electrodes on which we'll run ripple detection on, using RippleLFPSelection.set_lfp_electrodes

In [2]:

?sgr.RippleLFPSelection.set_lfp_electrodes

?sgr.RippleLFPSelection.set_lfp_electrodes

Signature:
sgrip.RippleLFPSelection.set_lfp_electrodes(
    key,
    electrode_list=None,
    group_name='CA1',
    **kwargs,
)
Docstring:
Removes all electrodes for the specified nwb file and then
adds back the electrodes in the list

Parameters
----------
key : dict
    dictionary corresponding to the LFPBand entry to use for
    ripple detection
electrode_list : list
    list of electrodes from LFPBandSelection.LFPBandElectrode
    to be used as the ripple LFP during detection
group_name : str, optional
    description of the electrode group, by default "CA1"
File:      ~/Src2/spyglass/src/spyglass/ripple/v1/ripple.py
Type:      function

We'll need the nwb_file_name, an electrode_list, and to a group_name.

• By default, group_name is set to CA1 for ripple detection, but we could alternatively use PFC.
• We use nwb_file_name to explore which electrodes are available for the electrode_list.

In [3]:

nwb_file_name = "tonks20211103_.nwb"
interval_list_name = "test interval"
filter_name = "Ripple 150-250 Hz"

nwb_file_name = "tonks20211103_.nwb" interval_list_name = "test interval" filter_name = "Ripple 150-250 Hz"

Now we can look at electrode_id in the Electrode table:

In [4]:

electrodes = (
    (sgc.Electrode() & {"nwb_file_name": nwb_file_name})
    * (
        lfp_analysis.LFPBandSelection.LFPBandElectrode()
        & {
            "nwb_file_name": nwb_file_name,
            "filter_name": filter_name,
            "target_interval_list_name": interval_list_name,
        }
    )
    * sgc.BrainRegion
).fetch(format="frame")
electrodes

electrodes = ( (sgc.Electrode() & {"nwb_file_name": nwb_file_name}) * ( lfp_analysis.LFPBandSelection.LFPBandElectrode() & { "nwb_file_name": nwb_file_name, "filter_name": filter_name, "target_interval_list_name": interval_list_name, } ) * sgc.BrainRegion ).fetch(format="frame") electrodes

Out[4]:

											probe_id	probe_shank	probe_electrode	name	original_reference_electrode	x	y	z	filtering	impedance	bad_channel	x_warped	y_warped	z_warped	contacts	region_name	subregion_name	subsubregion_name
nwb_file_name	electrode_group_name	electrode_id	lfp_merge_id	filter_name	filter_sampling_rate	target_interval_list_name	lfp_band_sampling_rate	lfp_electrode_group_name	reference_elect_id	region_id
tonks20211103_.nwb	7	28	2f3c93d5-5d5d-2d47-75b3-c346dddbd312	Ripple 150-250 Hz	1000	test interval	100	CA1_test	-1	19	tetrode_12.5	0	0	28	44	0.0	0.0	0.0	None	0.0	False	0.0	0.0	0.0		ca1	None	None
							1000	CA1_test	-1	19	tetrode_12.5	0	0	28	44	0.0	0.0	0.0	None	0.0	False	0.0	0.0	0.0		ca1	None	None
	8	32	2f3c93d5-5d5d-2d47-75b3-c346dddbd312	Ripple 150-250 Hz	1000	test interval	100	CA1_test	-1	19	tetrode_12.5	0	0	32	44	0.0	0.0	0.0	None	0.0	False	0.0	0.0	0.0		ca1	None	None
							1000	CA1_test	-1	19	tetrode_12.5	0	0	32	44	0.0	0.0	0.0	None	0.0	False	0.0	0.0	0.0		ca1	None	None

For ripple detection, we want only tetrodes, and only the first good wire on each tetrode. We will assume that is the first wire on each tetrode. I will do this using pandas syntax but you could use datajoint to filter this table as well. Here is the filtered table.

In [5]:

hpc_names = ["ca1", "hippocampus", "CA1", "Hippocampus"]
electrodes.loc[
    (electrodes.region_name.isin(hpc_names)) & (electrodes.probe_electrode == 0)
]

hpc_names = ["ca1", "hippocampus", "CA1", "Hippocampus"] electrodes.loc[ (electrodes.region_name.isin(hpc_names)) & (electrodes.probe_electrode == 0) ]

Out[5]:

											probe_id	probe_shank	probe_electrode	name	original_reference_electrode	x	y	z	filtering	impedance	bad_channel	x_warped	y_warped	z_warped	contacts	region_name	subregion_name	subsubregion_name
nwb_file_name	electrode_group_name	electrode_id	lfp_merge_id	filter_name	filter_sampling_rate	target_interval_list_name	lfp_band_sampling_rate	lfp_electrode_group_name	reference_elect_id	region_id
tonks20211103_.nwb	7	28	2f3c93d5-5d5d-2d47-75b3-c346dddbd312	Ripple 150-250 Hz	1000	test interval	100	CA1_test	-1	19	tetrode_12.5	0	0	28	44	0.0	0.0	0.0	None	0.0	False	0.0	0.0	0.0		ca1	None	None
							1000	CA1_test	-1	19	tetrode_12.5	0	0	28	44	0.0	0.0	0.0	None	0.0	False	0.0	0.0	0.0		ca1	None	None
	8	32	2f3c93d5-5d5d-2d47-75b3-c346dddbd312	Ripple 150-250 Hz	1000	test interval	100	CA1_test	-1	19	tetrode_12.5	0	0	32	44	0.0	0.0	0.0	None	0.0	False	0.0	0.0	0.0		ca1	None	None
							1000	CA1_test	-1	19	tetrode_12.5	0	0	32	44	0.0	0.0	0.0	None	0.0	False	0.0	0.0	0.0		ca1	None	None

We only want the electrode_id to put in the electrode_list:

In [13]:

electrode_list = np.unique(
    (
        electrodes.loc[
            (electrodes.region_name.isin(hpc_names))
            & (electrodes.probe_electrode == 0)
        ]
        .reset_index()
        .electrode_id
    ).tolist()
)

electrode_list.sort()

electrode_list = np.unique( ( electrodes.loc[ (electrodes.region_name.isin(hpc_names)) & (electrodes.probe_electrode == 0) ] .reset_index() .electrode_id ).tolist() ) electrode_list.sort()

By default, set_lfp_electrodes will use all the available electrodes from LFPBandV1.

We can insert into RippleLFPSelection and the RippleLFPElectrode part table, passing the key for the entry from LFPBandV1, our electrode_list, and the group_name into set_lfp_electrodes

In [18]:

group_name = "CA1_test"

lfp_band_key = (
    lfp_analysis.LFPBandV1()
    & {"filter_name": filter_name, "nwb_file_name": nwb_file_name}
).fetch1("KEY")

sgr.RippleLFPSelection.set_lfp_electrodes(
    lfp_band_key,
    electrode_list=electrode_list,
    group_name=group_name,
)

group_name = "CA1_test" lfp_band_key = ( lfp_analysis.LFPBandV1() & {"filter_name": filter_name, "nwb_file_name": nwb_file_name} ).fetch1("KEY") sgr.RippleLFPSelection.set_lfp_electrodes( lfp_band_key, electrode_list=electrode_list, group_name=group_name, )

In [20]:

sgr.RippleLFPSelection.RippleLFPElectrode()

sgr.RippleLFPSelection.RippleLFPElectrode()

Out[20]:

lfp_merge_id	filter_name descriptive name of this filter	filter_sampling_rate sampling rate for this filter	nwb_file_name name of the NWB file	target_interval_list_name descriptive name of this interval list	lfp_band_sampling_rate the sampling rate for this band	group_name	lfp_electrode_group_name the name of this group of electrodes	electrode_group_name electrode group name from NWBFile	electrode_id the unique number for this electrode	reference_elect_id the reference electrode to use; -1 for no reference
2f3c93d5-5d5d-2d47-75b3-c346dddbd312	Ripple 150-250 Hz	1000	tonks20211103_.nwb	test interval	100	CA1_test	CA1_test	7	28	-1
2f3c93d5-5d5d-2d47-75b3-c346dddbd312	Ripple 150-250 Hz	1000	tonks20211103_.nwb	test interval	100	CA1_test	CA1_test	8	32	-1

Total: 2

Here's the ripple selection key we'll use downstream

In [28]:

rip_sel_key = (sgrip.RippleLFPSelection & lfp_band_key).fetch1("KEY")

rip_sel_key = (sgrip.RippleLFPSelection & lfp_band_key).fetch1("KEY")

Setting Ripple Parameters¶

In [21]:

sgr.RippleParameters()

sgr.RippleParameters()

Out[21]:

ripple_param_name a name for this set of parameters	ripple_param_dict dictionary of parameters
default	=BLOB=

Total: 1

Here are the default ripple parameters:

In [22]:

(sgrip.RippleParameters() & {"ripple_param_name": "default"}).fetch1()

(sgrip.RippleParameters() & {"ripple_param_name": "default"}).fetch1()

Out[22]:

{'ripple_param_name': 'default',
 'ripple_param_dict': {'speed_name': 'head_speed',
  'ripple_detection_algorithm': 'Kay_ripple_detector',
  'ripple_detection_params': {'speed_threshold': 4.0,
   'minimum_duration': 0.015,
   'zscore_threshold': 2.0,
   'smoothing_sigma': 0.004,
   'close_ripple_threshold': 0.0}}}

• filter_name: which bandpass filter is used
• speed_name: the name of the speed parameters in IntervalPositionInfo

For the Kay_ripple_detector (options are currently Kay and Karlsson, see ripple_detection package for specifics) the parameters are:

• speed_threshold (cm/s): maximum speed the animal can move
• minimum_duration (s): minimum time above threshold
• zscore_threshold (std): minimum value to be considered a ripple, in standard deviations from mean
• smoothing_sigma (s): how much to smooth the signal in time
• close_ripple_threshold (s): exclude ripples closer than this amount

Check interval speed¶

The speed for this interval should exist under the default position parameter set and for a given interval.

In [27]:

pos_key = sgp.PositionOutput.merge_get_part(
    {
        "nwb_file_name": nwb_file_name,
        "position_info_param_name": "default",
        "interval_list_name": "pos 1 valid times",
    }
).fetch1("KEY")
(sgp.PositionOutput & pos_key).fetch1_dataframe()

pos_key = sgp.PositionOutput.merge_get_part( { "nwb_file_name": nwb_file_name, "position_info_param_name": "default", "interval_list_name": "pos 1 valid times", } ).fetch1("KEY") (sgp.PositionOutput & pos_key).fetch1_dataframe()

Out[27]:

	head_position_x	head_position_y	head_orientation	head_velocity_x	head_velocity_y	head_speed
time
1.635961e+09	98.670000	78.320000	1.878849	-0.212384	-1.050933e+00	1.072179
1.635961e+09	98.615000	78.210000	1.899349	-0.143244	-1.136351e+00	1.145344
1.635961e+09	98.633333	78.173333	1.919567	-0.031501	-1.123425e+00	1.123867
1.635961e+09	98.596667	78.100000	1.932884	0.094982	-1.013202e+00	1.017644
1.635961e+09	98.633333	78.100000	1.946067	0.194273	-8.272934e-01	0.849798
...	...	...	...	...	...	...
1.635963e+09	96.323333	71.500000	-2.265535	-0.415082	-1.486577e-05	0.415082
1.635963e+09	96.286667	71.500000	-2.158799	-0.413708	-3.243187e-06	0.413708
1.635963e+09	96.250000	71.500000	-2.034444	-0.374655	-6.383825e-07	0.374655
1.635963e+09	96.250000	71.500000	-2.034444	-0.307793	-1.133319e-07	0.307793
1.635963e+09	96.250000	71.500000	-2.034444	-0.229237	-1.813955e-08	0.229237

48950 rows × 6 columns

We'll use the head_speed above as part of RippleParameters.

Run Ripple Detection¶

Now we can put everything together.

In [32]:

key = {
    "ripple_param_name": "default",
    **rip_sel_key,
    "pos_merge_id": pos_key["merge_id"],
}
sgrip.RippleTimesV1().populate(key)

key = { "ripple_param_name": "default", **rip_sel_key, "pos_merge_id": pos_key["merge_id"], } sgrip.RippleTimesV1().populate(key)

Computing ripple times for: {'lfp_merge_id': UUID('2f3c93d5-5d5d-2d47-75b3-c346dddbd312'), 'filter_name': 'Ripple 150-250 Hz', 'filter_sampling_rate': 1000, 'nwb_file_name': 'tonks20211103_.nwb', 'target_interval_list_name': 'test interval', 'lfp_band_sampling_rate': 100, 'group_name': 'CA1_test', 'ripple_param_name': 'default', 'pos_merge_id': UUID('68959dc8-f8a3-c3c0-a534-096b3bc10f6c')}

<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
/home/dgramling/anaconda3/envs/spyglass-position3/lib/python3.9/site-packages/ripple_detection/detectors.py:31: RuntimeWarning: invalid value encountered in sqrt
  return np.sqrt(ripple_consensus_trace)
/home/dgramling/anaconda3/envs/spyglass-position3/lib/python3.9/site-packages/hdmf/spec/namespace.py:531: UserWarning: Ignoring cached namespace 'hdmf-common' version 1.5.1 because version 1.6.0 is already loaded.
  warn("Ignoring cached namespace '%s' version %s because version %s is already loaded."
/home/dgramling/anaconda3/envs/spyglass-position3/lib/python3.9/site-packages/hdmf/spec/namespace.py:531: UserWarning: Ignoring cached namespace 'core' version 2.4.0 because version 2.6.0-alpha is already loaded.
  warn("Ignoring cached namespace '%s' version %s because version %s is already loaded."
/home/dgramling/anaconda3/envs/spyglass-position3/lib/python3.9/site-packages/hdmf/spec/namespace.py:531: UserWarning: Ignoring cached namespace 'hdmf-experimental' version 0.2.0 because version 0.3.0 is already loaded.
  warn("Ignoring cached namespace '%s' version %s because version %s is already loaded."
/home/dgramling/anaconda3/envs/spyglass-position3/lib/python3.9/site-packages/pynwb/behavior.py:46: UserWarning: SpatialSeries 'series_0' has data shape (66792, 4) which is not compliant with NWB 2.5 and greater. The second dimension should have length <= 3 to represent at most x, y, z.
  warnings.warn("SpatialSeries '%s' has data shape %s which is not compliant with NWB 2.5 and greater. "
/home/dgramling/anaconda3/envs/spyglass-position3/lib/python3.9/site-packages/pynwb/behavior.py:46: UserWarning: SpatialSeries 'series_1' has data shape (48950, 4) which is not compliant with NWB 2.5 and greater. The second dimension should have length <= 3 to represent at most x, y, z.
  warnings.warn("SpatialSeries '%s' has data shape %s which is not compliant with NWB 2.5 and greater. "
/home/dgramling/anaconda3/envs/spyglass-position3/lib/python3.9/site-packages/pynwb/behavior.py:46: UserWarning: SpatialSeries 'series_2' has data shape (98507, 4) which is not compliant with NWB 2.5 and greater. The second dimension should have length <= 3 to represent at most x, y, z.
  warnings.warn("SpatialSeries '%s' has data shape %s which is not compliant with NWB 2.5 and greater. "
/home/dgramling/anaconda3/envs/spyglass-position3/lib/python3.9/site-packages/pynwb/behavior.py:46: UserWarning: SpatialSeries 'series_3' has data shape (44892, 4) which is not compliant with NWB 2.5 and greater. The second dimension should have length <= 3 to represent at most x, y, z.
  warnings.warn("SpatialSeries '%s' has data shape %s which is not compliant with NWB 2.5 and greater. "
/home/dgramling/anaconda3/envs/spyglass-position3/lib/python3.9/site-packages/pynwb/behavior.py:46: UserWarning: SpatialSeries 'series_4' has data shape (82313, 4) which is not compliant with NWB 2.5 and greater. The second dimension should have length <= 3 to represent at most x, y, z.
  warnings.warn("SpatialSeries '%s' has data shape %s which is not compliant with NWB 2.5 and greater. "
/home/dgramling/anaconda3/envs/spyglass-position3/lib/python3.9/site-packages/pynwb/behavior.py:46: UserWarning: SpatialSeries 'series_5' has data shape (81566, 4) which is not compliant with NWB 2.5 and greater. The second dimension should have length <= 3 to represent at most x, y, z.
  warnings.warn("SpatialSeries '%s' has data shape %s which is not compliant with NWB 2.5 and greater. "
/home/dgramling/anaconda3/envs/spyglass-position3/lib/python3.9/site-packages/pynwb/behavior.py:46: UserWarning: SpatialSeries 'series_6' has data shape (83811, 4) which is not compliant with NWB 2.5 and greater. The second dimension should have length <= 3 to represent at most x, y, z.
  warnings.warn("SpatialSeries '%s' has data shape %s which is not compliant with NWB 2.5 and greater. "

Writing new NWB file tonks20211103_6VF2MQ65RR.nwb

/home/dgramling/Src2/datajoint-python/datajoint/external.py:276: DeprecationWarning: The truth value of an empty array is ambiguous. Returning False, but in future this will result in an error. Use `array.size > 0` to check that an array is not empty.
  if check_hash:

And then fetch1_dataframe for ripple times

In [33]:

ripple_times = (sgrip.RippleTimesV1() & key).fetch1_dataframe()
ripple_times

ripple_times = (sgrip.RippleTimesV1() & key).fetch1_dataframe() ripple_times

Out[33]:

	start_time	end_time
id
0	1.635961e+09	1.635961e+09
1	1.635961e+09	1.635961e+09
2	1.635961e+09	1.635961e+09
3	1.635961e+09	1.635961e+09
4	1.635961e+09	1.635961e+09
5	1.635961e+09	1.635961e+09
6	1.635961e+09	1.635961e+09
7	1.635961e+09	1.635961e+09
8	1.635961e+09	1.635961e+09
9	1.635961e+09	1.635961e+09
10	1.635961e+09	1.635961e+09
11	1.635961e+09	1.635961e+09
12	1.635961e+09	1.635961e+09
13	1.635961e+09	1.635961e+09
14	1.635961e+09	1.635961e+09
15	1.635961e+09	1.635961e+09
16	1.635961e+09	1.635961e+09

Up Next¶

Next, we'll extract mark indicator.