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..a0f9202e --- /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.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/statistics.rst b/docs/api/statistics.rst index 1ab0ec89..b9761166 100644 --- a/docs/api/statistics.rst +++ b/docs/api/statistics.rst @@ -27,6 +27,15 @@ Spikestamps Statistics interspike_intervals coefficient_variation +Burst Analysis +------------------ + +.. autosummary:: + :nosignatures: + :toctree: _toctree/StatisticsAPI + + burst + Useful External Packages ======================== diff --git a/docs/api/visualization.rst b/docs/api/visualization.rst index d2836990..1edee119 100644 --- a/docs/api/visualization.rst +++ b/docs/api/visualization.rst @@ -46,6 +46,19 @@ Causality Analysis pairwise_causality_plot +Burst Analysis +------------------ + +.. currentmodule:: miv.visualization + +.. automodule:: miv.visualization.event + + .. autosummary:: + :nosignatures: + :toctree: _toctree/VisualizationAPI + + plot_burst + Useful External Packages ======================== diff --git a/docs/guide/burst_analysis.md b/docs/guide/burst_analysis.md new file mode 100644 index 00000000..9fbdcfb8 --- /dev/null +++ b/docs/guide/burst_analysis.md @@ -0,0 +1,80 @@ +--- +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][1],[2][2] + +[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 + ++++ + +## 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 +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 +burst(spiketrains,45,0.1,10) +``` + +## 3. Plotting +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/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 diff --git a/miv/statistics/__init__.py b/miv/statistics/__init__.py index ce34df9f..b07d0832 100644 --- a/miv/statistics/__init__.py +++ b/miv/statistics/__init__.py @@ -1,3 +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 new file mode 100644 index 00000000..bf6666b3 --- /dev/null +++ b/miv/statistics/burst.py @@ -0,0 +1,86 @@ +import os + +import matplotlib.pyplot as plt +import numpy as np + +from miv.statistics.spiketrain_statistics import interspike_intervals +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 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 = 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 + burst = [] # List to store burst parameters + + 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: + 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[:, 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) diff --git a/miv/visualization/__init__.py b/miv/visualization/__init__.py index 1678cdaf..755cdd6c 100644 --- a/miv/visualization/__init__.py +++ b/miv/visualization/__init__.py @@ -1,3 +1,4 @@ from miv.visualization.causality import * +from miv.visualization.event import * from miv.visualization.fft_domain import * from miv.visualization.waveform import * diff --git a/miv/visualization/event.py b/miv/visualization/event.py new file mode 100644 index 00000000..0bb872b4 --- /dev/null +++ b/miv/visualization/event.py @@ -0,0 +1,44 @@ +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 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]) 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