Burst analysis for a single channel and plotting across channels

docs/api/_toctree/StatisticsAPI/miv.statistics.burst.rst

@@ -0,0 +1,6 @@
+miv.statistics.burst
+====================
+
+.. currentmodule:: miv.statistics
+
+.. autofunction:: burst

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

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
 ========================
 

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
 ========================
 

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)
+
+```

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

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 *

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)

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 *

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

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])

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