From d4a703cd3170d0344f2ebb172db430fd6bbc11d6 Mon Sep 17 00:00:00 2001 From: Gauravu2 <89496329+Gauravu2@users.noreply.github.com> Date: Tue, 21 Jun 2022 10:01:12 -0500 Subject: [PATCH 01/17] Burst analysis for a single channel --- miv/statistics/__init__.py | 2 + miv/statistics/burst.py | 84 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 86 insertions(+) create mode 100644 miv/statistics/burst.py diff --git a/miv/statistics/__init__.py b/miv/statistics/__init__.py index ce34df9f..f5663f56 100644 --- a/miv/statistics/__init__.py +++ b/miv/statistics/__init__.py @@ -1,3 +1,5 @@ +from mib.statistics.burst import * + from miv.statistics.pairwise_causality import * from miv.statistics.signal_statistics import * from miv.statistics.spiketrain_statistics import * diff --git a/miv/statistics/burst.py b/miv/statistics/burst.py new file mode 100644 index 00000000..160911b5 --- /dev/null +++ b/miv/statistics/burst.py @@ -0,0 +1,84 @@ +# Bursting_Analysis: Burst length, location, and rate + +##Bursting is defined as the occurence of a specified number of spikes (usually >10), with a small interspike interval (usually < 100ms) + +import os + +import matplotlib.pyplot as plt +import numpy as np + +from miv.typing import SpikestampsType + + +def burst(spiketrains: SpikestampsType, channel: float, min_isi: float, min_len: float): + + # Calculates parameters critical to characterize bursting phenomenon on a singl channel + + # Parameters + # ---------- + # spikes : SpikestampsType + # Single spike-stamps + # Channel : float + # Channel to analyze + # min_isi : float + # Minimum Interspike Interval (in seconds) to be considered as bursting [standard = 0.1] + # min_len : float + # Minimum number of simultaneous spikes to be considered as bursting [standard = 10] + + # Returns + # ------- + # start_time: float + # The time instances when a burst starts + + # burst_duration: float + # The time duration of a particular burst + + # burst_len: float + # Number of spikes in a particular burst + + # burst_rate: float + # firing rates corresponding to particular bursts + + spike_interval = np.diff( + spiketrains[channel].magnitude + ) ## Calculate Inter Spike Interval (ISI) + A = np.array(spike_interval) + B = np.array(spike_interval) + B[A > min_isi] = 0 + B[A <= min_isi] = 1 ##Only spikes within specified min ISI are 1 otherwise 0 + PP = np.copy(B) + + Min_Spikes = min_len + P = [] + + for i in np.arange( + len(B) - Min_Spikes + ): ## Loop to check clusters of spikes greater than specified minimum length of burst + t = 0 + if np.sum(B[i : i + Min_Spikes]) == Min_Spikes: + q = 1 + while q > 0 and i + q + Min_Spikes <= len(B) - 1: + if B[i + q + Min_Spikes] == 1: + q += 1 + t = q + else: + q = 0 + P.append([i, t + Min_Spikes + 1]) + B[i : i + t + Min_Spikes] = 0 ## Zeroing counted spikes to avoid recounting + + Q = np.array(P) + if np.sum(Q) == 0: + start_time = 0 + end_time = 0 + burst_duration = 0 + burst_rate = 0 + burst_len = 0 + else: + spike = np.array(spiketrains[channel].magnitude) + start_time = spike[Q[:, 0]] + end_time = spike[Q[:, 0] + Q[:, 1]] + burst_len = Q[:, 1] + burst_duration = end_time - start_time + burst_rate = burst_len / (burst_duration) + + return start_time, burst_duration, burst_len, burst_rate From 29d662d359f0c7318bb1ed99ca4d2e6ecb96be10 Mon Sep 17 00:00:00 2001 From: Gauravu2 <89496329+Gauravu2@users.noreply.github.com> Date: Tue, 21 Jun 2022 10:10:47 -0500 Subject: [PATCH 02/17] Update __init__.py --- miv/statistics/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/miv/statistics/__init__.py b/miv/statistics/__init__.py index f5663f56..a78500bf 100644 --- a/miv/statistics/__init__.py +++ b/miv/statistics/__init__.py @@ -1,4 +1,4 @@ -from mib.statistics.burst import * +from miv.statistics.burst import * from miv.statistics.pairwise_causality import * from miv.statistics.signal_statistics import * From 7bc819930054044b8ca334ffde03bf9de34c296d Mon Sep 17 00:00:00 2001 From: Gauravu2 <89496329+Gauravu2@users.noreply.github.com> Date: Tue, 21 Jun 2022 21:04:40 -0500 Subject: [PATCH 03/17] Add files via upload --- miv/statistics/__init__.py | 1 - miv/statistics/burst.py | 2 +- miv/visualization/__init__.py | 4 ++- miv/visualization/plot_burst.py | 46 +++++++++++++++++++++++++++++++++ 4 files changed, 50 insertions(+), 3 deletions(-) create mode 100644 miv/visualization/plot_burst.py diff --git a/miv/statistics/__init__.py b/miv/statistics/__init__.py index a78500bf..b07d0832 100644 --- a/miv/statistics/__init__.py +++ b/miv/statistics/__init__.py @@ -1,5 +1,4 @@ from miv.statistics.burst import * - from miv.statistics.pairwise_causality import * from miv.statistics.signal_statistics import * from miv.statistics.spiketrain_statistics import * diff --git a/miv/statistics/burst.py b/miv/statistics/burst.py index 160911b5..6c824c6d 100644 --- a/miv/statistics/burst.py +++ b/miv/statistics/burst.py @@ -76,7 +76,7 @@ def burst(spiketrains: SpikestampsType, channel: float, min_isi: float, min_len: else: spike = np.array(spiketrains[channel].magnitude) start_time = spike[Q[:, 0]] - end_time = spike[Q[:, 0] + Q[:, 1]] + end_time = spike[Q[:, 0] + Q[:, 1] - 1] burst_len = Q[:, 1] burst_duration = end_time - start_time burst_rate = burst_len / (burst_duration) diff --git a/miv/visualization/__init__.py b/miv/visualization/__init__.py index 1678cdaf..33b53725 100644 --- a/miv/visualization/__init__.py +++ b/miv/visualization/__init__.py @@ -1,3 +1,5 @@ -from miv.visualization.causality import * from miv.visualization.fft_domain import * +from miv.visualization.pairwise_causality_plot import * +from miv.visualization.plot_burst import * +from miv.visualization.plot_spectral import * from miv.visualization.waveform import * diff --git a/miv/visualization/plot_burst.py b/miv/visualization/plot_burst.py new file mode 100644 index 00000000..e20a1b0e --- /dev/null +++ b/miv/visualization/plot_burst.py @@ -0,0 +1,46 @@ +# Bursting_Plot: Plots bursts as a bar raster with corresponding start time and duration +##Bursting is defined as the occurence of a specified number of spikes (usually >10), with a small interspike interval (usually < 100ms) + + +import os + +import matplotlib.pyplot as plt +import numpy as np +from miv.statistics import burst +from miv.typing import SpikestampsType + + +def plot_burst(spiketrains: SpikestampsType, min_isi: float, min_len: float): + + # Plots burst events across electrodes to characterize bursting phenomenon on a singl channel + + # Parameters + # ---------- + # spikes : SpikestampsType + # Single spike-stamps + # min_isi : float + # Minimum Interspike Interval (in seconds) to be considered as bursting [standard = 0.1] + # min_len : float + # Minimum number of simultaneous spikes to be considered as bursting [standard = 10] + + # Returns + # ------- + # figure, axes + # matplot figure with bursts plotted for all electrodes + + fig, ax = plt.subplots() + start_time = [] + burst_duration = [] + burst_len = [] + burst_rate = [] + for i in np.arange(len(spiketrains)): + start_time, burst_duration, burst_len, burst_rate = burst( + spiketrains, i, min_isi, min_len + ) + a = np.column_stack((start_time, burst_duration)) + ax.broken_barh(a, (i + 1, 0.5), facecolors="tab:orange") + + ax.set_xlabel("Time (s)") + ax.set_ylabel("Electrode") + + return fig, ax From 2a97e5ff9ef554746d9a327d1bc850f887611713 Mon Sep 17 00:00:00 2001 From: Gauravu2 <89496329+Gauravu2@users.noreply.github.com> Date: Tue, 21 Jun 2022 21:06:31 -0500 Subject: [PATCH 04/17] Documentation for burst analysis and plotting --- docs/guide/burst_analysis.md | 171 +++++++++++++++++++++++++++++++++++ 1 file changed, 171 insertions(+) create mode 100644 docs/guide/burst_analysis.md diff --git a/docs/guide/burst_analysis.md b/docs/guide/burst_analysis.md new file mode 100644 index 00000000..8f8ae4de --- /dev/null +++ b/docs/guide/burst_analysis.md @@ -0,0 +1,171 @@ +--- +jupytext: + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.13.8 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# Burst Analysis + +Bursting is defined as the occurence of a specified number of simultaneuos spikes (usually >10), with a small interspike interval (usually < 100ms) [1,2] + +References: + +[1] Chiappalone, Michela, et al. "Burst detection algorithms for the analysis of spatio-temporal patterns in cortical networks of neurons." Neurocomputing 65 (2005): 653-662.\ + +[2] Eisenman, Lawrence N., et al. "Quantification of bursting and synchrony in cultured hippocampal neurons." Journal of neurophysiology 114.2 (2015): 1059-1071. + ++++ + +## 1. Data Load and Preprocessing + +```{code-cell} ipython3 +:tags: [hide-cell] + +import os, sys +import miv +import scipy.signal as ss + +from miv.io import DataManager +from miv.signal.filter import ButterBandpass, MedianFilter, FilterCollection +from miv.signal.spike import ThresholdCutoff +from miv.statistics import pairwise_causality +from miv.visualization import plot_spectral, pairwise_causality_plot +import numpy as np +from miv.typing import SignalType +``` + +```{code-cell} ipython3 +# Experiment name +experiment_query = "experiment0" + +# Data call +signal_filter = ( + FilterCollection() + .append(ButterBandpass(600, 2400, order=4)) + .append(MedianFilter(threshold=60, k=30)) +) +spike_detection = ThresholdCutoff(cutoff=5) + +# Spike Detection +data_collection = optogenetic.load_data() +data = data_collection.query_path_name(experiment_query)[0] +#false_channels = [12,15,36,41,42,43,45,47,53,57,55,58,61,62] +#data.set_channel_mask(false_channels) +with data.load() as (signal, timestamps, sampling_rate): + # Preprocess + signal = signal_filter(signal, sampling_rate) + spiketrains = spike_detection(signal, timestamps, sampling_rate) +``` + +## 2. Burst Estimations + + Function Call + + burst(spiketrains: SpikestampsType, channel: float, min_isi: float, min_len: float) + + Calculates parameters critical to characterize bursting phenomenon on a single channel + + Parameters + ---------- + spikes : SpikestampsType + Single spike-stamps + Channel : float + Channel to analyze + min_isi : float + Minimum Interspike Interval (in seconds) to be considered as bursting [standard = 0.1] + min_len : float + Minimum number of simultaneous spikes to be considered as bursting [standard = 10] + + Returns + ------- + start_time: float + The time instances when a burst starts + + burst_duration: float + The time duration of a particular burst + + burst_len: float + Number of spikes in a particular burst + + burst_rate: float + firing rates corresponding to particular bursts + + + +```{code-cell} ipython3 +#Example +burst(spiketrains,45,0.1,10) +# Estimates the burst parameters for 45th electrode with bursts defined as more than 10 simultaneous spikes with 0.1 s interspike interval +``` + +## 3. Plot bursting events across electrodes + + Function Call: + + plot_burst(spiketrains: SpikestampsType, min_isi: float, min_len: float) + + Parameters + ---------- + spikes : SpikestampsType + Single spike-stamps + min_isi : float + Minimum Interspike Interval (in seconds) to be considered as bursting [standard = 0.1] + min_len : float + Minimum number of simultaneous spikes to be considered as bursting [standard = 10] + + Returns + ------- + figure, axes + matplot figure with bursts plotted for all electrodes + + +```{code-cell} ipython3 +#Example +plot_burst(spiketrains,0.1,10) +# plots the burst events fo with bursts defined as more than 10 simultaneous spikes with 0.1 s interspike interval +``` + +## Welch Coherence + +Plots Power Spectral Densities for channels X and Y, Cross Power Spctral Densities and Coherence between them using Welch's method\ + +plot_spectral(signal, X, Y, sampling_rate, Number_Segments)\ + + +Parameters\ +#---------- +signal : SignalType\ + Input signal\ +X : float\ + First Channel \ +Y : float\ + Second Channel\ +sampling_rate : float\ + Sampling frequency\ +Number_Segments: float\ +Number of segments to divide the entire signal\ + +Returns\ +#------- +figure: plt.Figure\ +axes\ + +References:\ + +1) https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.coherence.html +2) P. Welch, “The use of the fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms”, IEEE Trans. Audio Electroacoust. vol. 15, pp. 70-73, 1967. + + +```{code-cell} ipython3 +#Example +plot_spectral(signal,1,42,30000,10000) + +##Plots the PSDs, CPSD and Coherence for channel 1 & 43 for a sampling rate of 30000, with signal divided into 10000 segments +``` From 469a0e26f35e4293c90c44998b035aa9e4e48ff8 Mon Sep 17 00:00:00 2001 From: Gauravu2 <89496329+Gauravu2@users.noreply.github.com> Date: Tue, 21 Jun 2022 21:17:54 -0500 Subject: [PATCH 05/17] updated init --- miv/visualization/__init__.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/miv/visualization/__init__.py b/miv/visualization/__init__.py index 33b53725..91870333 100644 --- a/miv/visualization/__init__.py +++ b/miv/visualization/__init__.py @@ -1,5 +1,4 @@ +from miv.visualization.causality import * from miv.visualization.fft_domain import * -from miv.visualization.pairwise_causality_plot import * -from miv.visualization.plot_burst import * -from miv.visualization.plot_spectral import * from miv.visualization.waveform import * +from miv.visualization.plot_burst import * \ No newline at end of file From d08c0927312bfdbdd00ec8bd3b11d194ed483652 Mon Sep 17 00:00:00 2001 From: Gauravu2 <89496329+Gauravu2@users.noreply.github.com> Date: Tue, 21 Jun 2022 21:50:34 -0500 Subject: [PATCH 06/17] Update burst analysis --- docs/guide/burst_analysis.md | 101 +++----------------------------- miv/statistics/burst.py | 68 ++++++++++----------- miv/visualization/plot_burst.py | 39 ++++++------ 3 files changed, 59 insertions(+), 149 deletions(-) diff --git a/docs/guide/burst_analysis.md b/docs/guide/burst_analysis.md index 8f8ae4de..bbfd8b78 100644 --- a/docs/guide/burst_analysis.md +++ b/docs/guide/burst_analysis.md @@ -64,108 +64,21 @@ with data.load() as (signal, timestamps, sampling_rate): spiketrains = spike_detection(signal, timestamps, sampling_rate) ``` -## 2. Burst Estimations - - Function Call - - burst(spiketrains: SpikestampsType, channel: float, min_isi: float, min_len: float) - - Calculates parameters critical to characterize bursting phenomenon on a single channel - - Parameters - ---------- - spikes : SpikestampsType - Single spike-stamps - Channel : float - Channel to analyze - min_isi : float - Minimum Interspike Interval (in seconds) to be considered as bursting [standard = 0.1] - min_len : float - Minimum number of simultaneous spikes to be considered as bursting [standard = 10] - - Returns - ------- - start_time: float - The time instances when a burst starts - - burst_duration: float - The time duration of a particular burst - - burst_len: float - Number of spikes in a particular burst - - burst_rate: float - firing rates corresponding to particular bursts - - +## 2. Burst Estimations +Calculates parameters critical to characterize bursting phenomenon on a single channel. Documentation is available [here](miv.statistics.burst). ```{code-cell} ipython3 -#Example -burst(spiketrains,45,0.1,10) # Estimates the burst parameters for 45th electrode with bursts defined as more than 10 simultaneous spikes with 0.1 s interspike interval +burst(spiketrains,45,0.1,10) ``` -## 3. Plot bursting events across electrodes - - Function Call: - - plot_burst(spiketrains: SpikestampsType, min_isi: float, min_len: float) - - Parameters - ---------- - spikes : SpikestampsType - Single spike-stamps - min_isi : float - Minimum Interspike Interval (in seconds) to be considered as bursting [standard = 0.1] - min_len : float - Minimum number of simultaneous spikes to be considered as bursting [standard = 10] - - Returns - ------- - figure, axes - matplot figure with bursts plotted for all electrodes - +## 3. Plotting +Plots the burst events across the recordings. Documentation is available [here](miv.visualization.plot_burst). ```{code-cell} ipython3 #Example +# plots the burst events with bursts defined as more than 10 simultaneous spikes with 0.1 s interspike interval plot_burst(spiketrains,0.1,10) -# plots the burst events fo with bursts defined as more than 10 simultaneous spikes with 0.1 s interspike interval + ``` -## Welch Coherence - -Plots Power Spectral Densities for channels X and Y, Cross Power Spctral Densities and Coherence between them using Welch's method\ - -plot_spectral(signal, X, Y, sampling_rate, Number_Segments)\ - - -Parameters\ -#---------- -signal : SignalType\ - Input signal\ -X : float\ - First Channel \ -Y : float\ - Second Channel\ -sampling_rate : float\ - Sampling frequency\ -Number_Segments: float\ -Number of segments to divide the entire signal\ - -Returns\ -#------- -figure: plt.Figure\ -axes\ - -References:\ - -1) https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.coherence.html -2) P. Welch, “The use of the fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms”, IEEE Trans. Audio Electroacoust. vol. 15, pp. 70-73, 1967. - - -```{code-cell} ipython3 -#Example -plot_spectral(signal,1,42,30000,10000) - -##Plots the PSDs, CPSD and Coherence for channel 1 & 43 for a sampling rate of 30000, with signal divided into 10000 segments -``` diff --git a/miv/statistics/burst.py b/miv/statistics/burst.py index 6c824c6d..01593722 100644 --- a/miv/statistics/burst.py +++ b/miv/statistics/burst.py @@ -1,7 +1,3 @@ -# Bursting_Analysis: Burst length, location, and rate - -##Bursting is defined as the occurence of a specified number of spikes (usually >10), with a small interspike interval (usually < 100ms) - import os import matplotlib.pyplot as plt @@ -11,34 +7,38 @@ def burst(spiketrains: SpikestampsType, channel: float, min_isi: float, min_len: float): - - # Calculates parameters critical to characterize bursting phenomenon on a singl channel - - # Parameters - # ---------- - # spikes : SpikestampsType - # Single spike-stamps - # Channel : float - # Channel to analyze - # min_isi : float - # Minimum Interspike Interval (in seconds) to be considered as bursting [standard = 0.1] - # min_len : float - # Minimum number of simultaneous spikes to be considered as bursting [standard = 10] - - # Returns - # ------- - # start_time: float - # The time instances when a burst starts - - # burst_duration: float - # The time duration of a particular burst - - # burst_len: float - # Number of spikes in a particular burst - - # burst_rate: float - # firing rates corresponding to particular bursts - + """ + Calculates parameters critical to characterize bursting phenomenon on a single channel + Bursting is defined as the occurence of a specified number of spikes (usually >10), with a small interspike interval (usually < 100ms) [1]_, [2]_ + + Parameters + ---------- + spikes : SpikestampsType + Single spike-stamps + Channel : float + Channel to analyze + min_isi : float + Minimum Interspike Interval (in seconds) to be considered as bursting [standard = 0.1] + min_len : float + Minimum number of simultaneous spikes to be considered as bursting [standard = 10] + + Returns + ------- + start_time: float + The time instances when a burst starts + burst_duration: float + The time duration of a particular burst + burst_len: float + Number of spikes in a particular burst + burst_rate: float + firing rates corresponding to particular bursts + + ..[1] Chiappalone, Michela, et al. "Burst detection algorithms for the analysis of spatio-temporal patterns + in cortical networks of neurons." Neurocomputing 65 (2005): 653-662. + ..[2] Eisenman, Lawrence N., et al. "Quantification of bursting and synchrony in cultured + hippocampal neurons." Journal of neurophysiology 114.2 (2015): 1059-1071. + + """ spike_interval = np.diff( spiketrains[channel].magnitude ) ## Calculate Inter Spike Interval (ISI) @@ -63,7 +63,7 @@ def burst(spiketrains: SpikestampsType, channel: float, min_isi: float, min_len: t = q else: q = 0 - P.append([i, t + Min_Spikes + 1]) + P.append([i, t + Min_Spikes]) B[i : i + t + Min_Spikes] = 0 ## Zeroing counted spikes to avoid recounting Q = np.array(P) @@ -76,7 +76,7 @@ def burst(spiketrains: SpikestampsType, channel: float, min_isi: float, min_len: else: spike = np.array(spiketrains[channel].magnitude) start_time = spike[Q[:, 0]] - end_time = spike[Q[:, 0] + Q[:, 1] - 1] + end_time = spike[Q[:, 0] + Q[:, 1]] burst_len = Q[:, 1] burst_duration = end_time - start_time burst_rate = burst_len / (burst_duration) diff --git a/miv/visualization/plot_burst.py b/miv/visualization/plot_burst.py index e20a1b0e..c7aa6edd 100644 --- a/miv/visualization/plot_burst.py +++ b/miv/visualization/plot_burst.py @@ -1,7 +1,3 @@ -# Bursting_Plot: Plots bursts as a bar raster with corresponding start time and duration -##Bursting is defined as the occurence of a specified number of spikes (usually >10), with a small interspike interval (usually < 100ms) - - import os import matplotlib.pyplot as plt @@ -11,23 +7,24 @@ def plot_burst(spiketrains: SpikestampsType, min_isi: float, min_len: float): - - # Plots burst events across electrodes to characterize bursting phenomenon on a singl channel - - # Parameters - # ---------- - # spikes : SpikestampsType - # Single spike-stamps - # min_isi : float - # Minimum Interspike Interval (in seconds) to be considered as bursting [standard = 0.1] - # min_len : float - # Minimum number of simultaneous spikes to be considered as bursting [standard = 10] - - # Returns - # ------- - # figure, axes - # matplot figure with bursts plotted for all electrodes - + + """ + Plots burst events across electrodes to characterize bursting phenomenon on a singl channel + + Parameters + ---------- + spikes : SpikestampsType + Single spike-stamps + min_isi : float + Minimum Interspike Interval (in seconds) to be considered as bursting [standard = 0.1] + min_len : float + Minimum number of simultaneous spikes to be considered as bursting [standard = 10] + + Returns + ------- + figure, axes + matplot figure with bursts plotted for all electrodes + """ fig, ax = plt.subplots() start_time = [] burst_duration = [] From 748ed277c9c5e93d97f15638727687a3ff5ff04e Mon Sep 17 00:00:00 2001 From: Gauravu2 <89496329+Gauravu2@users.noreply.github.com> Date: Wed, 22 Jun 2022 22:50:20 -0500 Subject: [PATCH 07/17] doc: update burst_analysis --- docs/guide/burst_analysis.md | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/docs/guide/burst_analysis.md b/docs/guide/burst_analysis.md index bbfd8b78..c9430eef 100644 --- a/docs/guide/burst_analysis.md +++ b/docs/guide/burst_analysis.md @@ -13,13 +13,10 @@ kernelspec: # Burst Analysis -Bursting is defined as the occurence of a specified number of simultaneuos spikes (usually >10), with a small interspike interval (usually < 100ms) [1,2] +Bursting is defined as the occurence of a specified number of simultaneuos spikes (usually >10), with a small interspike interval (usually < 100ms) [1][1],[2][2] -References: - -[1] Chiappalone, Michela, et al. "Burst detection algorithms for the analysis of spatio-temporal patterns in cortical networks of neurons." Neurocomputing 65 (2005): 653-662.\ - -[2] Eisenman, Lawrence N., et al. "Quantification of bursting and synchrony in cultured hippocampal neurons." Journal of neurophysiology 114.2 (2015): 1059-1071. +[1]: https://www.sciencedirect.com/science/article/abs/pii/S0925231204004874 +[2]: https://journals.physiology.org/doi/full/10.1152/jn.00079.2015?rfr_dat=cr_pub++0pubmed&url_ver=Z39.88-2003&rfr_id=ori%3Arid%3Acrossref.org +++ From b1aec22c4ab3ef04ffd5ff7ba528e1790122da37 Mon Sep 17 00:00:00 2001 From: Gauravu2 <89496329+Gauravu2@users.noreply.github.com> Date: Wed, 22 Jun 2022 23:07:32 -0500 Subject: [PATCH 08/17] update: algorithm burst.py --- miv/statistics/burst.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/miv/statistics/burst.py b/miv/statistics/burst.py index 01593722..b1c60a57 100644 --- a/miv/statistics/burst.py +++ b/miv/statistics/burst.py @@ -7,7 +7,7 @@ def burst(spiketrains: SpikestampsType, channel: float, min_isi: float, min_len: float): - """ + """ Calculates parameters critical to characterize bursting phenomenon on a single channel Bursting is defined as the occurence of a specified number of spikes (usually >10), with a small interspike interval (usually < 100ms) [1]_, [2]_ @@ -38,10 +38,8 @@ def burst(spiketrains: SpikestampsType, channel: float, min_isi: float, min_len: ..[2] Eisenman, Lawrence N., et al. "Quantification of bursting and synchrony in cultured hippocampal neurons." Journal of neurophysiology 114.2 (2015): 1059-1071. - """ - spike_interval = np.diff( - spiketrains[channel].magnitude - ) ## Calculate Inter Spike Interval (ISI) + """ + spike_interval = np.diff(spiketrains[channel].magnitude) ## Calculate Inter Spike Interval (ISI) A = np.array(spike_interval) B = np.array(spike_interval) B[A > min_isi] = 0 From d130d3bc87f2a66d8e9fd1b4306e41eae65e97cd Mon Sep 17 00:00:00 2001 From: Gauravu2 <89496329+Gauravu2@users.noreply.github.com> Date: Wed, 22 Jun 2022 23:16:01 -0500 Subject: [PATCH 09/17] Update plot_burst.py --- miv/visualization/plot_burst.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/miv/visualization/plot_burst.py b/miv/visualization/plot_burst.py index c7aa6edd..1816b84b 100644 --- a/miv/visualization/plot_burst.py +++ b/miv/visualization/plot_burst.py @@ -8,7 +8,7 @@ def plot_burst(spiketrains: SpikestampsType, min_isi: float, min_len: float): - """ + """ Plots burst events across electrodes to characterize bursting phenomenon on a singl channel Parameters From e5ab06fccf424566d0c3485bdc356c49a2b09e90 Mon Sep 17 00:00:00 2001 From: Gauravu2 <89496329+Gauravu2@users.noreply.github.com> Date: Wed, 22 Jun 2022 23:34:05 -0500 Subject: [PATCH 10/17] doc: update for burst and plotting usage --- .../_toctree/StatisticsAPI/miv.statistics.burst.rst | 6 ++++++ .../miv.visualization.plot_burst.rst | 6 ++++++ docs/api/statistics.rst | 13 +++++++++++++ docs/api/visualization.rst | 13 +++++++++++++ 4 files changed, 38 insertions(+) create mode 100644 docs/api/_toctree/StatisticsAPI/miv.statistics.burst.rst create mode 100644 docs/api/_toctree/VisualizationAPI/miv.visualization.plot_burst.rst diff --git a/docs/api/_toctree/StatisticsAPI/miv.statistics.burst.rst b/docs/api/_toctree/StatisticsAPI/miv.statistics.burst.rst new file mode 100644 index 00000000..1941cb65 --- /dev/null +++ b/docs/api/_toctree/StatisticsAPI/miv.statistics.burst.rst @@ -0,0 +1,6 @@ + miv.statistics.burst + ===================================================== + + .. currentmodule:: miv.statistics + + .. autofunction:: burst diff --git a/docs/api/_toctree/VisualizationAPI/miv.visualization.plot_burst.rst b/docs/api/_toctree/VisualizationAPI/miv.visualization.plot_burst.rst new file mode 100644 index 00000000..6f2deda7 --- /dev/null +++ b/docs/api/_toctree/VisualizationAPI/miv.visualization.plot_burst.rst @@ -0,0 +1,6 @@ + miv.visualization.plot\_burst + ===================================================== + + .. currentmodule:: miv.visualization + + .. autofunction:: plot_burst diff --git a/docs/api/statistics.rst b/docs/api/statistics.rst index 1ab0ec89..32b57fd3 100644 --- a/docs/api/statistics.rst +++ b/docs/api/statistics.rst @@ -26,6 +26,19 @@ Spikestamps Statistics firing_rates interspike_intervals coefficient_variation + +Burst Analysis +------------------ + +.. currentmodule:: miv.statistics + +.. automodule:: miv.statistics.burst + +.. autosummary:: + :nosignatures: + :toctree: _toctree/StatisticsAPI + + burst Useful External Packages ======================== diff --git a/docs/api/visualization.rst b/docs/api/visualization.rst index d2836990..20e54dfd 100644 --- a/docs/api/visualization.rst +++ b/docs/api/visualization.rst @@ -45,6 +45,19 @@ Causality Analysis :toctree: _toctree/VisualizationAPI pairwise_causality_plot + +Burst Analysis +------------------ + +.. currentmodule:: miv.visualization + +.. automodule:: miv.visualization.plot_burst + + .. autosummary:: + :nosignatures: + :toctree: _toctree/VisualizationAPI + + plot_burst Useful External Packages ======================== From 94ce1d98c8d97799fca885d78114a40d798bb807 Mon Sep 17 00:00:00 2001 From: Gauravu2 <89496329+Gauravu2@users.noreply.github.com> Date: Thu, 23 Jun 2022 01:22:06 -0500 Subject: [PATCH 11/17] formatting burst and plotting --- miv/statistics/burst.py | 62 +++++++++++++++++---------------- miv/visualization/__init__.py | 2 +- miv/visualization/plot_burst.py | 33 +++++++++--------- 3 files changed, 50 insertions(+), 47 deletions(-) diff --git a/miv/statistics/burst.py b/miv/statistics/burst.py index b1c60a57..7f4a15a7 100644 --- a/miv/statistics/burst.py +++ b/miv/statistics/burst.py @@ -8,38 +8,40 @@ def burst(spiketrains: SpikestampsType, channel: float, min_isi: float, min_len: float): """ - Calculates parameters critical to characterize bursting phenomenon on a single channel - Bursting is defined as the occurence of a specified number of spikes (usually >10), with a small interspike interval (usually < 100ms) [1]_, [2]_ - - Parameters - ---------- - spikes : SpikestampsType - Single spike-stamps - Channel : float - Channel to analyze - min_isi : float - Minimum Interspike Interval (in seconds) to be considered as bursting [standard = 0.1] - min_len : float - Minimum number of simultaneous spikes to be considered as bursting [standard = 10] + Calculates parameters critical to characterize bursting phenomenon on a single channel + Bursting is defined as the occurence of a specified number of spikes (usually >10), with a small interspike interval (usually < 100ms) [1]_, [2]_ + + Parameters + ---------- + spikes : SpikestampsType + Single spike-stamps + Channel : float + Channel to analyze + min_isi : float + Minimum Interspike Interval (in seconds) to be considered as bursting [standard = 0.1] + min_len : float + Minimum number of simultaneous spikes to be considered as bursting [standard = 10] + + Returns + ------- + start_time: float + The time instances when a burst starts + burst_duration: float + The time duration of a particular burst + burst_len: float + Number of spikes in a particular burst + burst_rate: float + firing rates corresponding to particular bursts + + ..[1] Chiappalone, Michela, et al. "Burst detection algorithms for the analysis of spatio-temporal patterns + in cortical networks of neurons." Neurocomputing 65 (2005): 653-662. + ..[2] Eisenman, Lawrence N., et al. "Quantification of bursting and synchrony in cultured + hippocampal neurons." Journal of neurophysiology 114.2 (2015): 1059-1071. - Returns - ------- - start_time: float - The time instances when a burst starts - burst_duration: float - The time duration of a particular burst - burst_len: float - Number of spikes in a particular burst - burst_rate: float - firing rates corresponding to particular bursts - - ..[1] Chiappalone, Michela, et al. "Burst detection algorithms for the analysis of spatio-temporal patterns - in cortical networks of neurons." Neurocomputing 65 (2005): 653-662. - ..[2] Eisenman, Lawrence N., et al. "Quantification of bursting and synchrony in cultured - hippocampal neurons." Journal of neurophysiology 114.2 (2015): 1059-1071. - """ - spike_interval = np.diff(spiketrains[channel].magnitude) ## Calculate Inter Spike Interval (ISI) + spike_interval = np.diff( + spiketrains[channel].magnitude + ) ## Calculate Inter Spike Interval (ISI) A = np.array(spike_interval) B = np.array(spike_interval) B[A > min_isi] = 0 diff --git a/miv/visualization/__init__.py b/miv/visualization/__init__.py index 91870333..45712431 100644 --- a/miv/visualization/__init__.py +++ b/miv/visualization/__init__.py @@ -1,4 +1,4 @@ from miv.visualization.causality import * from miv.visualization.fft_domain import * +from miv.visualization.plot_burst import * from miv.visualization.waveform import * -from miv.visualization.plot_burst import * \ No newline at end of file diff --git a/miv/visualization/plot_burst.py b/miv/visualization/plot_burst.py index 1816b84b..9113079a 100644 --- a/miv/visualization/plot_burst.py +++ b/miv/visualization/plot_burst.py @@ -2,28 +2,29 @@ import matplotlib.pyplot as plt import numpy as np + from miv.statistics import burst from miv.typing import SpikestampsType def plot_burst(spiketrains: SpikestampsType, min_isi: float, min_len: float): - + """ - Plots burst events across electrodes to characterize bursting phenomenon on a singl channel - - Parameters - ---------- - spikes : SpikestampsType - Single spike-stamps - min_isi : float - Minimum Interspike Interval (in seconds) to be considered as bursting [standard = 0.1] - min_len : float - Minimum number of simultaneous spikes to be considered as bursting [standard = 10] - - Returns - ------- - figure, axes - matplot figure with bursts plotted for all electrodes + Plots burst events across electrodes to characterize bursting phenomenon on a singl channel + + Parameters + ---------- + spikes : SpikestampsType + Single spike-stamps + min_isi : float + Minimum Interspike Interval (in seconds) to be considered as bursting [standard = 0.1] + min_len : float + Minimum number of simultaneous spikes to be considered as bursting [standard = 10] + + Returns + ------- + figure, axes + matplot figure with bursts plotted for all electrodes """ fig, ax = plt.subplots() start_time = [] From 65721d67f084ff34478913f9c1e2f8c846a2cb95 Mon Sep 17 00:00:00 2001 From: Gauravu2 <89496329+Gauravu2@users.noreply.github.com> Date: Thu, 23 Jun 2022 14:16:25 -0500 Subject: [PATCH 12/17] test: add test case for burst --- miv/statistics/burst.py | 39 +++++++++++++++++----------------- tests/statistics/test_burst.py | 27 +++++++++++++++++++++++ 2 files changed, 47 insertions(+), 19 deletions(-) create mode 100644 tests/statistics/test_burst.py diff --git a/miv/statistics/burst.py b/miv/statistics/burst.py index 7f4a15a7..8469205f 100644 --- a/miv/statistics/burst.py +++ b/miv/statistics/burst.py @@ -39,34 +39,35 @@ def burst(spiketrains: SpikestampsType, channel: float, min_isi: float, min_len: hippocampal neurons." Journal of neurophysiology 114.2 (2015): 1059-1071. """ + spike_interval = np.diff( spiketrains[channel].magnitude - ) ## Calculate Inter Spike Interval (ISI) - A = np.array(spike_interval) - B = np.array(spike_interval) - B[A > min_isi] = 0 - B[A <= min_isi] = 1 ##Only spikes within specified min ISI are 1 otherwise 0 - PP = np.copy(B) - - Min_Spikes = min_len - P = [] + ) # Calculate Inter Spike Interval (ISI) + burst_spike = (spike_interval <= min_isi).astype( + np.int32 + ) # Only spikes within specified min ISI are 1 otherwise 0 and are stored + min_spikes = min_len + burst = [] # List to store burst parameters for i in np.arange( - len(B) - Min_Spikes - ): ## Loop to check clusters of spikes greater than specified minimum length of burst + len(burst_spike) - min_spikes + ): # Loop to check clusters of spikes greater than specified minimum length of burst t = 0 - if np.sum(B[i : i + Min_Spikes]) == Min_Spikes: + if np.sum(burst_spike[i : i + min_spikes]) == min_spikes: q = 1 - while q > 0 and i + q + Min_Spikes <= len(B) - 1: - if B[i + q + Min_Spikes] == 1: - q += 1 + while q > 0 and i + q + min_spikes <= len(burst_spike) - 1: + if burst_spike[i + q + min_spikes - 1] == 1: t = q + q += 1 else: q = 0 - P.append([i, t + Min_Spikes]) - B[i : i + t + Min_Spikes] = 0 ## Zeroing counted spikes to avoid recounting - Q = np.array(P) + burst.append([i, t + min_spikes]) + burst_spike[ + i : i + t + min_spikes + ] = 0 # Zeroing counted spikes to avoid recounting + Q = np.array(burst) + if np.sum(Q) == 0: start_time = 0 end_time = 0 @@ -77,7 +78,7 @@ def burst(spiketrains: SpikestampsType, channel: float, min_isi: float, min_len: spike = np.array(spiketrains[channel].magnitude) start_time = spike[Q[:, 0]] end_time = spike[Q[:, 0] + Q[:, 1]] - burst_len = Q[:, 1] + burst_len = Q[:, 1] + 1 burst_duration = end_time - start_time burst_rate = burst_len / (burst_duration) diff --git a/tests/statistics/test_burst.py b/tests/statistics/test_burst.py new file mode 100644 index 00000000..0a9a469a --- /dev/null +++ b/tests/statistics/test_burst.py @@ -0,0 +1,27 @@ +import numpy as np +import pytest +from neo.core import AnalogSignal, Segment, SpikeTrain + +# Test set For Burst Module + + +def test_burst_analysis_output(): + from miv.statistics import burst + + # Initialize the spiketrain as below + seg = Segment(index=1) + train0 = SpikeTrain( + times=[0.1, 1.2, 1.3, 1.4, 1.5, 1.6, 4, 5, 5.1, 5.2, 8, 9.5], + units="sec", + t_stop=10, + ) + # train0 = SpikeTrain(times=[0.1,4,5,5.1,5.2,9.5], units='sec', t_stop=10) + seg.spiketrains.append(train0) + + output = burst(seg.spiketrains, 0, 0.2, 2) + # The function above should return two burst events from 1.2 (duration 0.4, length 5, rate 12.5 ) + # and 5(duration 0.2, length 3, rate 15) + np.testing.assert_allclose(output[0], [1.2, 5]) + np.testing.assert_allclose(output[1], [0.4, 0.2]) + np.testing.assert_allclose(output[2], [5, 3]) + np.testing.assert_allclose(output[3], [12.5, 15]) From d57dd7565f062381c2ab022defc2fc187048a7ef Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Fri, 24 Jun 2022 04:29:58 -0500 Subject: [PATCH 13/17] update: add burst doc --- docs/index.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/index.rst b/docs/index.rst index 9b7deb63..a9f02e52 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -39,6 +39,7 @@ If you are interested in contributing to this project, we prepared contribution guide/spike_sorting guide/channel_correlation guide/granger_causality_psd_cpsd_coherence + guide/burst_analysis .. toctree:: :maxdepth: 2 From 3d146231182f03b6e84427832ba2deb56dd130ab Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Fri, 24 Jun 2022 04:31:03 -0500 Subject: [PATCH 14/17] fix: documentation build errors --- docs/api/_toctree/StatisticsAPI/miv.statistics.burst.rst | 8 ++++---- .../miv.visualization.event.plot_burst.rst | 6 ++++++ .../VisualizationAPI/miv.visualization.plot_burst.rst | 6 ------ docs/api/statistics.rst | 6 +----- docs/api/visualization.rst | 4 ++-- docs/guide/burst_analysis.md | 9 ++++----- miv/statistics/burst.py | 4 ++-- miv/visualization/__init__.py | 2 +- miv/visualization/{plot_burst.py => event.py} | 2 +- 9 files changed, 21 insertions(+), 26 deletions(-) create mode 100644 docs/api/_toctree/VisualizationAPI/miv.visualization.event.plot_burst.rst delete mode 100644 docs/api/_toctree/VisualizationAPI/miv.visualization.plot_burst.rst rename miv/visualization/{plot_burst.py => event.py} (100%) diff --git a/docs/api/_toctree/StatisticsAPI/miv.statistics.burst.rst b/docs/api/_toctree/StatisticsAPI/miv.statistics.burst.rst index 1941cb65..a0f9202e 100644 --- a/docs/api/_toctree/StatisticsAPI/miv.statistics.burst.rst +++ b/docs/api/_toctree/StatisticsAPI/miv.statistics.burst.rst @@ -1,6 +1,6 @@ - miv.statistics.burst - ===================================================== +miv.statistics.burst +==================== - .. currentmodule:: miv.statistics +.. currentmodule:: miv.statistics - .. autofunction:: burst +.. autofunction:: burst diff --git a/docs/api/_toctree/VisualizationAPI/miv.visualization.event.plot_burst.rst b/docs/api/_toctree/VisualizationAPI/miv.visualization.event.plot_burst.rst new file mode 100644 index 00000000..f019b4e4 --- /dev/null +++ b/docs/api/_toctree/VisualizationAPI/miv.visualization.event.plot_burst.rst @@ -0,0 +1,6 @@ +miv.visualization.event.plot\_burst +=================================== + +.. currentmodule:: miv.visualization.event + +.. autofunction:: plot_burst diff --git a/docs/api/_toctree/VisualizationAPI/miv.visualization.plot_burst.rst b/docs/api/_toctree/VisualizationAPI/miv.visualization.plot_burst.rst deleted file mode 100644 index 6f2deda7..00000000 --- a/docs/api/_toctree/VisualizationAPI/miv.visualization.plot_burst.rst +++ /dev/null @@ -1,6 +0,0 @@ - miv.visualization.plot\_burst - ===================================================== - - .. currentmodule:: miv.visualization - - .. autofunction:: plot_burst diff --git a/docs/api/statistics.rst b/docs/api/statistics.rst index 32b57fd3..b9761166 100644 --- a/docs/api/statistics.rst +++ b/docs/api/statistics.rst @@ -26,14 +26,10 @@ Spikestamps Statistics firing_rates interspike_intervals coefficient_variation - + Burst Analysis ------------------ -.. currentmodule:: miv.statistics - -.. automodule:: miv.statistics.burst - .. autosummary:: :nosignatures: :toctree: _toctree/StatisticsAPI diff --git a/docs/api/visualization.rst b/docs/api/visualization.rst index 20e54dfd..1edee119 100644 --- a/docs/api/visualization.rst +++ b/docs/api/visualization.rst @@ -45,13 +45,13 @@ Causality Analysis :toctree: _toctree/VisualizationAPI pairwise_causality_plot - + Burst Analysis ------------------ .. currentmodule:: miv.visualization -.. automodule:: miv.visualization.plot_burst +.. automodule:: miv.visualization.event .. autosummary:: :nosignatures: diff --git a/docs/guide/burst_analysis.md b/docs/guide/burst_analysis.md index c9430eef..9fbdcfb8 100644 --- a/docs/guide/burst_analysis.md +++ b/docs/guide/burst_analysis.md @@ -61,21 +61,20 @@ with data.load() as (signal, timestamps, sampling_rate): spiketrains = spike_detection(signal, timestamps, sampling_rate) ``` -## 2. Burst Estimations +## 2. Burst Estimations Calculates parameters critical to characterize bursting phenomenon on a single channel. Documentation is available [here](miv.statistics.burst). ```{code-cell} ipython3 -# Estimates the burst parameters for 45th electrode with bursts defined as more than 10 simultaneous spikes with 0.1 s interspike interval +# Estimates the burst parameters for 45th electrode with bursts defined as more than 10 simultaneous spikes with 0.1 s interspike interval burst(spiketrains,45,0.1,10) ``` ## 3. Plotting -Plots the burst events across the recordings. Documentation is available [here](miv.visualization.plot_burst). +Plots the burst events across the recordings. Documentation is available [here](miv.visualization.event.plot_burst). ```{code-cell} ipython3 #Example # plots the burst events with bursts defined as more than 10 simultaneous spikes with 0.1 s interspike interval plot_burst(spiketrains,0.1,10) - -``` +``` diff --git a/miv/statistics/burst.py b/miv/statistics/burst.py index 8469205f..b046ff6e 100644 --- a/miv/statistics/burst.py +++ b/miv/statistics/burst.py @@ -33,9 +33,9 @@ def burst(spiketrains: SpikestampsType, channel: float, min_isi: float, min_len: burst_rate: float firing rates corresponding to particular bursts - ..[1] Chiappalone, Michela, et al. "Burst detection algorithms for the analysis of spatio-temporal patterns + .. [1] Chiappalone, Michela, et al. "Burst detection algorithms for the analysis of spatio-temporal patterns in cortical networks of neurons." Neurocomputing 65 (2005): 653-662. - ..[2] Eisenman, Lawrence N., et al. "Quantification of bursting and synchrony in cultured + .. [2] Eisenman, Lawrence N., et al. "Quantification of bursting and synchrony in cultured hippocampal neurons." Journal of neurophysiology 114.2 (2015): 1059-1071. """ diff --git a/miv/visualization/__init__.py b/miv/visualization/__init__.py index 45712431..755cdd6c 100644 --- a/miv/visualization/__init__.py +++ b/miv/visualization/__init__.py @@ -1,4 +1,4 @@ from miv.visualization.causality import * +from miv.visualization.event import * from miv.visualization.fft_domain import * -from miv.visualization.plot_burst import * from miv.visualization.waveform import * diff --git a/miv/visualization/plot_burst.py b/miv/visualization/event.py similarity index 100% rename from miv/visualization/plot_burst.py rename to miv/visualization/event.py index 9113079a..0bb872b4 100644 --- a/miv/visualization/plot_burst.py +++ b/miv/visualization/event.py @@ -8,7 +8,6 @@ def plot_burst(spiketrains: SpikestampsType, min_isi: float, min_len: float): - """ Plots burst events across electrodes to characterize bursting phenomenon on a singl channel @@ -26,6 +25,7 @@ def plot_burst(spiketrains: SpikestampsType, min_isi: float, min_len: float): figure, axes matplot figure with bursts plotted for all electrodes """ + fig, ax = plt.subplots() start_time = [] burst_duration = [] From 80d4368a51fb8dc0942d66a8ade7e31f442ea34a Mon Sep 17 00:00:00 2001 From: Seung Hyun Kim Date: Sat, 25 Jun 2022 23:54:24 -0500 Subject: [PATCH 15/17] optimize burst algorithm --- miv/statistics/burst.py | 61 +++++++++++++++++++++-------------------- 1 file changed, 32 insertions(+), 29 deletions(-) diff --git a/miv/statistics/burst.py b/miv/statistics/burst.py index b046ff6e..a89d86df 100644 --- a/miv/statistics/burst.py +++ b/miv/statistics/burst.py @@ -3,6 +3,7 @@ import matplotlib.pyplot as plt import numpy as np +from miv.statistics.spiketrain_statistics import interspike_intervals from miv.typing import SpikestampsType @@ -14,7 +15,7 @@ def burst(spiketrains: SpikestampsType, channel: float, min_isi: float, min_len: Parameters ---------- spikes : SpikestampsType - Single spike-stamps + Single spike-stamps Channel : float Channel to analyze min_isi : float @@ -25,13 +26,13 @@ def burst(spiketrains: SpikestampsType, channel: float, min_isi: float, min_len: Returns ------- start_time: float - The time instances when a burst starts + The time instances when a burst starts burst_duration: float - The time duration of a particular burst + The time duration of a particular burst burst_len: float - Number of spikes in a particular burst + Number of spikes in a particular burst burst_rate: float - firing rates corresponding to particular bursts + firing rates corresponding to particular bursts .. [1] Chiappalone, Michela, et al. "Burst detection algorithms for the analysis of spatio-temporal patterns in cortical networks of neurons." Neurocomputing 65 (2005): 653-662. @@ -40,32 +41,23 @@ def burst(spiketrains: SpikestampsType, channel: float, min_isi: float, min_len: """ - spike_interval = np.diff( - spiketrains[channel].magnitude - ) # Calculate Inter Spike Interval (ISI) + spike_interval = interspike_intervals(spiketrains[channel].magnitude) burst_spike = (spike_interval <= min_isi).astype( - np.int32 + np.bool_ ) # Only spikes within specified min ISI are 1 otherwise 0 and are stored - min_spikes = min_len burst = [] # List to store burst parameters - for i in np.arange( - len(burst_spike) - min_spikes - ): # Loop to check clusters of spikes greater than specified minimum length of burst - t = 0 - if np.sum(burst_spike[i : i + min_spikes]) == min_spikes: - q = 1 - while q > 0 and i + q + min_spikes <= len(burst_spike) - 1: - if burst_spike[i + q + min_spikes - 1] == 1: - t = q - q += 1 - else: - q = 0 - - burst.append([i, t + min_spikes]) - burst_spike[ - i : i + t + min_spikes - ] = 0 # Zeroing counted spikes to avoid recounting + flag = False + start_idx = -1 + delta = np.logical_xor(burst_spike[:-1], burst_spike[1:]) + for idx, dval in enumerate(delta): + if dval: + flag = ~flag + if flag: + start_idx = idx + 1 + else: + if idx + 1 - start_idx >= min_len: + burst.append((start_idx, idx + 1)) Q = np.array(burst) if np.sum(Q) == 0: @@ -77,9 +69,20 @@ def burst(spiketrains: SpikestampsType, channel: float, min_isi: float, min_len: else: spike = np.array(spiketrains[channel].magnitude) start_time = spike[Q[:, 0]] - end_time = spike[Q[:, 0] + Q[:, 1]] - burst_len = Q[:, 1] + 1 + end_time = spike[Q[:, 1]] + burst_len = Q[:, 1] - Q[:, 0] + 1 burst_duration = end_time - start_time burst_rate = burst_len / (burst_duration) return start_time, burst_duration, burst_len, burst_rate + + +if __name__ == "__main__": + import timeit + + from neo.core import SpikeTrain + + arr = np.random.random(10000) + arr = np.sort(arr) + train0 = [SpikeTrain(times=arr, units="sec", t_stop=arr.max())] + o_algorithm = timeit.timeit(lambda: burst(train0, 0, 0.2, 5), number=1000) From f5ff3721af5fa49d72247a4c79d62b6760bc6cf6 Mon Sep 17 00:00:00 2001 From: Gauravu2 <89496329+Gauravu2@users.noreply.github.com> Date: Fri, 24 Jun 2022 12:34:56 -0500 Subject: [PATCH 16/17] Updated to assert error --- miv/statistics/burst.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/miv/statistics/burst.py b/miv/statistics/burst.py index a89d86df..bf6666b3 100644 --- a/miv/statistics/burst.py +++ b/miv/statistics/burst.py @@ -11,7 +11,6 @@ def burst(spiketrains: SpikestampsType, channel: float, min_isi: float, min_len: """ Calculates parameters critical to characterize bursting phenomenon on a single channel Bursting is defined as the occurence of a specified number of spikes (usually >10), with a small interspike interval (usually < 100ms) [1]_, [2]_ - Parameters ---------- spikes : SpikestampsType @@ -22,7 +21,6 @@ def burst(spiketrains: SpikestampsType, channel: float, min_isi: float, min_len: Minimum Interspike Interval (in seconds) to be considered as bursting [standard = 0.1] min_len : float Minimum number of simultaneous spikes to be considered as bursting [standard = 10] - Returns ------- start_time: float @@ -38,10 +36,10 @@ def burst(spiketrains: SpikestampsType, channel: float, min_isi: float, min_len: in cortical networks of neurons." Neurocomputing 65 (2005): 653-662. .. [2] Eisenman, Lawrence N., et al. "Quantification of bursting and synchrony in cultured hippocampal neurons." Journal of neurophysiology 114.2 (2015): 1059-1071. - """ spike_interval = interspike_intervals(spiketrains[channel].magnitude) + assert spike_interval.all() > 0, "Inter Spike Interval cannot be zero" burst_spike = (spike_interval <= min_isi).astype( np.bool_ ) # Only spikes within specified min ISI are 1 otherwise 0 and are stored From f60495f985dcad80998c346b7386da59b231c933 Mon Sep 17 00:00:00 2001 From: Gauravu2 <89496329+Gauravu2@users.noreply.github.com> Date: Fri, 24 Jun 2022 12:35:31 -0500 Subject: [PATCH 17/17] Updated test case --- tests/test_burst.py | 38 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) create mode 100644 tests/test_burst.py diff --git a/tests/test_burst.py b/tests/test_burst.py new file mode 100644 index 00000000..61428f1b --- /dev/null +++ b/tests/test_burst.py @@ -0,0 +1,38 @@ +import numpy as np +import pytest +from neo.core import AnalogSignal, Segment, SpikeTrain + +# Test set For Burst Module + + +def test_burst_analysis_output(): + from miv.statistics import burst + + # Initialize the spiketrain as below + seg = Segment(index=1) + train0 = SpikeTrain( + times=[0.1, 1.2, 1.3, 1.4, 1.5, 1.6, 4, 5, 5.1, 5.2, 8, 9.5], + units="sec", + t_stop=10, + ) + # train0 = SpikeTrain(times=[0.1,4,5,5.1,5.2,9.5], units='sec', t_stop=10) + seg.spiketrains.append(train0) + + output = burst(seg.spiketrains, 0, 0.2, 2) + # The function above should return two burst events from 1.2 (duration 0.4, length 5, rate 12.5 ) + # and 5(duration 0.2, length 3, rate 15) + np.testing.assert_allclose(output[0], [1.2, 5]) + np.testing.assert_allclose(output[1], [0.4, 0.2]) + np.testing.assert_allclose(output[2], [5, 3]) + np.testing.assert_allclose(output[3], [12.5, 15]) + + seg1 = Segment(index=1) + train1 = SpikeTrain( + times=[0.1, 1.2, 1.2, 1.2, 5, 6], + units="sec", + t_stop=10, + ) + seg1.spiketrains.append(train1) + with np.testing.assert_raises(AssertionError): + output = burst(seg1.spiketrains, 0, 0.2, 2) + # The function above should throw an error since the rate will become 1/0