# ~/analyseur/cbgtc/stat/sync.py
#
# Documentation by Lungsi 16 Oct 2025
#
import numpy as np
from analyseur.cbgtc.curate import get_desired_spiketimes_subset
from analyseur.cbgtc.parameters import SignalAnalysisParams
from analyseur.cbgtc.stats.rate import Rate
from analyseur.cbgtc.stats.psth import PSTH
from scipy.stats import pearsonr
[docs]
class Synchrony(object):
"""
Computes measures of synchrony among the neurons with given spike times
+----------------------------------+----------------------------------------------------------------------------------------------------------+
| Methods | Argument |
+==================================+==========================================================================================================+
| :py:meth:`.compute_basic` | - `spiketimes_superset`: see :class:`~analyseur.cbgtc.loader.LoadSpikeTimes.get_spiketimes_superset` |
| | - OPTIONAL: `binsz` (0.01 [default]), `window` ((0, 10) [default]) |
+----------------------------------+----------------------------------------------------------------------------------------------------------+
| :py:meth:`.compute_basic_slide` | - `spiketimes_superset`: see :class:`~analyseur.cbgtc.loader.LoadSpikeTimes.get_spiketimes_superset` |
| | - OPTIONAL: `binsz` (0.01 [default]), `window` ((0, 10) [default]), `windowsz` (0.5 [default]) |
+----------------------------------+----------------------------------------------------------------------------------------------------------+
| :py:meth:`.compute_fano_factor` | - `spiketimes_superset`: see :class:`~analyseur.cbgtc.loader.LoadSpikeTimes.get_spiketimes_superset` |
| | - OPTIONAL: `binsz` (0.01 [default]), `window` ((0, 10) [default]) |
+----------------------------------+----------------------------------------------------------------------------------------------------------+
=========
Use Cases
=========
------------------
1. Pre-requisites
------------------
1.1. Import Modules
````````````````````
::
from analyseur.cbgtc.loader import LoadSpikeTimes
from analyseur.cbgtc.stats.sync import Synchrony
1.2. Load file and get spike times
```````````````````````````````````
::
loadST = LoadSpikeTimes("spikes_GPi.csv")
spiketimes_superset = loadST.get_spiketimes_superset()
---------
2. Cases
---------
2.1. Compute basic measure of spike times synchrony (for all neurons)
``````````````````````````````````````````````````````````````````````
::
B = Synchrony.compute_basic(spiketimes_superset)
This returns the value for
.. math::
Sync = \\sqrt{\\frac{var\\left(\\left[\\mu\\left(\\left[f^{{i}}(t)\\right]_{\\forall{t}}\\right)\\right]_{\\forall{i}}\\right)}{\\mu\\left(\\left[var\\left(\\left[f^{(i)}(t)\\right]_{\\forall{i}}\\right)\\right]_{\\forall{t}}\\right)}}
2.2. Compute the basic measure of synchrony on a smoother frequency estimation
```````````````````````````````````````````````````````````````````````````````
::
S = Synchrony.compute_basic_slide(spiketimes_superset)
This returns the value for
.. math::
Sync = \\sqrt{\\frac{var\\left(\\left[\\mu\\left(\\left[f^{{i}}(t)\\right]_{\\forall{t}}\\right)\\right]_{\\forall{i}}\\right)}{\\mu\\left(\\left[var\\left(\\left[f^{(i)}(t)\\right]_{\\forall{i}}\\right)\\right]_{\\forall{t}}\\right)}}
2.3. Compute Fano factor as a metric for measuring spike times synchrony (for all neurons)
```````````````````````````````````````````````````````````````````````````````````````````
::
Fs = Synchrony.compute_fano_factor(spiketimes_superset)
This returns the value for
.. math::
F_{Sync} = \\frac{var\\left(\\left[\\sum_{\\forall{i}}p^{(i)}(t)\\right]_{\\forall{t}}\\right)}{\\mu\\left(\\left[\\sum_{\\forall{i}}p^{(i)}(t)\\right]_{\\forall{t}}\\right)}
.. raw:: html
<hr style="border: 2px solid red; margin: 20px 0;">
"""
__siganal = SignalAnalysisParams()
@staticmethod
def __get_spike_matrix(spiketimes_set, window, binsz):
[desired_spiketimes_subset, _] = get_desired_spiketimes_subset(spiketimes_set, neurons="all")
n_neurons = len(desired_spiketimes_subset)
time_bins = np.arange(window[0], window[1] + binsz, binsz)
n_bins = len(time_bins)
spike_matrix = np.zeros((n_neurons, n_bins))
# Fill the spike matrix
for i, indiv_spiketimes in enumerate(desired_spiketimes_subset):
for spiketime in indiv_spiketimes:
if window[0] <= spiketime < window[1]:
j = int((spiketime - window[0]) // binsz)
spike_matrix[i, j] = 1
return spike_matrix, time_bins
@staticmethod
def __count_matrix_and_active_neurons(spiketimes_set, window, binsz):
min_spikes = 1
n_neurons = len(spiketimes_set)
duration = window[1] - window[0]
# Create time bins
n_bins = int(duration / binsz)
t_bins = np.linspace(0, duration, n_bins)
# Create spike count matrix and identify active neurons
spike_counts = np.zeros((n_neurons, n_bins-1))
active_neurons = []
for i, (nX, indiv_spiketimes) in enumerate(spiketimes_set.items()):
counts = np.histogram(indiv_spiketimes, bins=t_bins)[0]
if np.sum(counts) >= min_spikes and np.std(counts) > 0:
spike_counts[i, :] = counts
active_neurons.append(i)
return spike_counts, active_neurons, t_bins
@staticmethod
def __get_spikearray_and_window(spiketimes_superset, window, neurons="all"):
[desired_spiketimes_subset, _] = get_desired_spiketimes_subset(spiketimes_superset, neurons=neurons)
n_neurons = len(desired_spiketimes_subset)
spike_arrays = [np.array(spike_times) for spike_times in desired_spiketimes_subset]
if window is None:
all_spikes = np.concatenate(spike_arrays)
window = (np.min(all_spikes), np.max(all_spikes))
return spike_arrays, window
@staticmethod
def __compute_for_basic(freq_matrix, filter_zeros=False):
# if filter_zeros: # population frequencies during active periods, i.e, excluding time bins with no activity
# freq_matrix = freq_matrix[freq_matrix > 0]
# freq_matrix[:, np.all(freq_matrix > 0, axis=0)]
if filter_zeros:
n = freq_matrix.shape[1]
colmn_wise_means = np.zeros(n)
colmn_wise_vars = np.zeros(n)
for j in range(n):
freqs = freq_matrix[:,j]
# Mean of rates across neurons for all t
colmn_wise_means[j] = np.mean(freqs[freqs > 0])
# Variance of rates across neurons for all t
colmn_wise_vars[j] = np.var(freqs[freqs > 0])
else:
# Mean of rates across t
colmn_wise_means = np.mean(freq_matrix, axis=0)
# Variance of rates across t
colmn_wise_vars = np.var(freq_matrix, axis=0)
variance_F = np.var(colmn_wise_means)
mean_F = np.mean(colmn_wise_vars)
if mean_F == 0:
if variance_F == 0:
s_sync = 0.0 # Edge Case: Perfect synchrony
else:
s_sync = np.inf # Theoretical: Infinite synchrony
else:
s_sync = np.sqrt(variance_F / mean_F)
return s_sync, colmn_wise_means, colmn_wise_vars
@staticmethod
def __compute_fano(pop_spike_count_matrix, filter_zeros=False):
# if filter_spikes: # population spike counts during active periods, i.e, excluding time bins with no activity
# pop_spike_count_matrix = pop_spike_count_matrix[pop_spike_count_matrix > 0]
# if len(pop_spike_count_matrix) > 0:
# # Sum of spike counts across neurons for all t
# colmn_wise_sums = np.sum(pop_spike_count_matrix, axis=0)
#
# variance_S = np.var(colmn_wise_sums)
# mean_S = np.mean(colmn_wise_sums)
#
# if mean_S == 0:
# fanofactor = 0.0
# else:
# fanofactor = variance_S / mean_S
# else:
# colmn_wise_sums = np.array([0])
# fanofactor = 0.0 # variance_S = 0.0, mean_S = 0.0
colmn_wise_sums = np.sum(pop_spike_count_matrix, axis=0)
if filter_zeros:
colmn_wise_sums = colmn_wise_sums[colmn_wise_sums > 0]
variance_S = np.var(colmn_wise_sums)
mean_S = np.mean(colmn_wise_sums)
if mean_S == 0:
fanofactor = 0.0
else:
fanofactor = variance_S / mean_S
return fanofactor, colmn_wise_sums
# @staticmethod
# def __compute_pairwise_corr(pop_spike_count_matrix):
# correlations = []
# pairs = []
#
# n_neurons = pop_spike_count_matrix.shape[0]
#
# for i in range(n_neurons):
# for j in range(i+1, n_neurons):
# corr, p_val = pearsonr(pop_spike_count_matrix[i], pop_spike_count_matrix[j])
# correlations.append(corr)
# pairs.append((i,j))
#
# return np.array(correlations), pairs
@staticmethod
def __compute_pairwise_corr(spike_count_matrix, active_neurons):
correlations = []
pairs = []
active_counts = spike_count_matrix[active_neurons, :]
n_active = len(active_neurons)
# for i in range(n_active):
# for j in range(i+1, n_active):
# # Check if either neuron has constant activity
# if np.std(active_counts[i]) > 0 and np.std(active_counts[j]) > 0:
# try:
# corr, p_val = pearsonr(active_counts[i], active_counts[j])
# correlations.append(corr)
# pairs.append((active_neurons[i], active_neurons[j]))
# except:
# correlations.append(0)
# pairs.append((0,0))
# Remove redundant std computations
stds = np.std(active_counts, axis=1) # compute once
for i in range(n_active):
for j in range(i+1, n_active):
if stds[i] > 0 and stds[j] > 0:
corr, p_val = pearsonr(active_counts[i], active_counts[j])
correlations.append(corr)
return np.array(correlations), pairs
[docs]
@classmethod
def compute_basic(cls, spiketimes_set, binsz=None, window=None, filter_zeros=True):
"""
Returns the basic measure of synchrony of spiking from all neurons.
:param spiketimes_set: Dictionary returned using :meth:`~analyseur.cbgtc.loader.LoadSpikeTimes.get_spiketimes_superset`
or using :meth:`~analyseur.cbgtc.loader.LoadSpikeTimes.get_spiketimes_subset`
:param binsz: integer or float; `0.01` [default]
:param window: Tuple in the form `(start_time, end_time)`; `(0, 10)` [default]
:param filter_zeros: boolean; `True` [default]
:return:
* s_sync : float
* freq_matrix : numpy array, Population frequencies per time bin
* time_bins_center : numpy array, Center times of each bin
**Formula**
.. table:: Formula
====================================================================================================================== ===============================================================
Definitions Interpretation
====================================================================================================================== ===============================================================
total neurons, :math:`n_{Nuc}` total number of neurons in the Nucleus
neuron index, :math:`i` i-th neuron in the pool of :math:`n_{Nuc}` neurons
frequency, :math:`f^{(i)}(t)` frequency of the i-th neuron at time :math:`t`
frequency matrix, :math:`F = \\left[f(a,b) = f^{(a)}(b)\\right]_{\\forall{a \\in [1, n_{Nuc}], b \\in [t_0, t_T]}}` frequencies of all (:math:`n_{Nuc}`) neurons for all times
====================================================================================================================== ===============================================================
Let the :math:`var(\\cdot)`, `variance function <https://numpy.org/doc/stable/reference/generated/numpy.var.html>`_ and
the :math:`\\mu(\\cdot)`, `arithmetic mean function <https://numpy.org/doc/stable/reference/generated/numpy.mean.html>`_
be implemented as shown
.. math::
F = \\overset{\\begin{matrix}t_0 & \\quad\\quad & t_1 & & & &\\ldots & & & t_T\\end{matrix}}
{\\underset{
\\begin{matrix}
\\quad\\quad\\uparrow & \\quad\\quad\\quad & \\uparrow & \\quad &\\ldots & & & \\uparrow \n
\\quad\\mu_{t_0} & \\quad\\quad\\quad & \\mu_{t_1} & \\quad &\\ldots & & & \\mu_{t_T} & \\rightarrow var_{\\forall{t}}
\\end{matrix}}
{\\begin{bmatrix}
f^{(1)}(t_0) & f^{(1)}(t_1) & \\ldots & f^{(1)}(t_T) \n
f^{(2)}(t_0) & f^{(2)}(t_1) & \\ldots & f^{(2)}(t_T) \n
\\vdots & \\vdots & \\ldots & \\vdots \n
f^{(i)}(t_0) & f^{(i)}(t_1) & \\ldots & f^{(i)}(t_T) \n
\\vdots & \\vdots & \\ldots & \\vdots \n
f^{(n_{Nuc})}(t_0) & f^{(n_{Nuc})}(t_1) & \\ldots & f^{(n_{Nuc})}(t_T)
\\end{bmatrix}
}}
Then, we define
.. math::
A \\triangleq var\\left(\\begin{bmatrix}
\\mu_{t_0} & \\mu_{t_1} & \\ldots & \\mu_{t_T}
\\end{bmatrix}\\right) = var_{\\forall{t}}
Implementing the `variance function <https://numpy.org/doc/stable/reference/generated/numpy.var.html>`_ and
the `arithmetic mean function <https://numpy.org/doc/stable/reference/generated/numpy.mean.html>`_ as shown below
.. math::
F = \\overset{\\begin{matrix}t_0 & \\quad\\quad & t_1 & & & &\\ldots & & & t_T\\end{matrix}}
{\\underset{
\\begin{matrix}
\\quad\\quad\\uparrow & \\quad\\quad\\quad & \\uparrow & \\quad &\\ldots & & & \\uparrow \n
\\quad var_{t_0} & \\quad\\quad\\quad & var_{t_1} & \\quad &\\ldots & & & var_{t_T} & \\rightarrow \\mu_{\\forall{t}}
\\end{matrix}}
{\\begin{bmatrix}
f^{(1)}(t_0) & f^{(1)}(t_1) & \\ldots & f^{(1)}(t_T) \n
f^{(2)}(t_0) & f^{(2)}(t_1) & \\ldots & f^{(2)}(t_T) \n
\\vdots & \\vdots & \\ldots & \\vdots \n
f^{(i)}(t_0) & f^{(i)}(t_1) & \\ldots & f^{(i)}(t_T) \n
\\vdots & \\vdots & \\ldots & \\vdots \n
f^{(n_{Nuc})}(t_0) & f^{(n_{Nuc})}(t_1) & \\ldots & f^{(n_{Nuc})}(t_T)
\\end{bmatrix}
}}
we make another definition
.. math::
B \\triangleq \\mu\\left(\\begin{bmatrix}
var_{t_0} & var_{t_1} & \\ldots & var_{t_T}
\\end{bmatrix}\\right) = \\mu_{\\forall{t}}
Then, synchrony is measured as
.. math::
Sync = \\sqrt{\\frac{A}{B}} = \\sqrt{\\frac{var\\left(\\left[\\mu\\left(\\left[f^{{i}}(t)\\right]_{\\forall{i}}\\right)\\right]_{\\forall{t}}\\right)}{\\mu\\left(\\left[var\\left(\\left[f^{(i)}(t)\\right]_{\\forall{i}}\\right)\\right]_{\\forall{t}}\\right)}}
NOTE: This method is a simple histogram-based approach that uses fixed bins.
**INTERPRETATION**
.. table:: Interpretation
======================= =============================================
:math:`S_{sync}` Interpretation
======================= =============================================
:math:`\\approx 0` Low synchrony (independent firing)
:math:`\\ge 1` High synchrony (neurons fire together)
:math:`= \\infty` Perfect synchrony
======================= =============================================
COMMENTS:
* Fano factor can indicate population-level synchrony
- higher values correspond to more synchrony
* Fano factor depends on time bin size.
**Notes on Filtering**
The default argument `filter_zeros=True` means only frequencies during active periods are accounted for
making the computation.
.. raw:: html
<hr style="border: 2px solid red; margin: 20px 0;">
"""
#============== DEFAULT Parameters ==============
if window is None:
window = cls.__siganal.window
if binsz is None:
binsz = cls.__siganal.binsz_100perbin
_, freq_matrix, time_bins = Rate.get_count_rate_matrix(spiketimes_set=spiketimes_set,
window=window, binsz=binsz,
neurons="all")
time_bins_center = (time_bins[:-1] + time_bins[1:]) / 2
[s_sync, _, _] = cls.__compute_for_basic(freq_matrix, filter_zeros=filter_zeros)
return s_sync, freq_matrix, time_bins_center
[docs]
@classmethod
def compute_basic_slide(cls, spiketimes_set, binsz=None, window=None, windowsz=None):
"""
Returns the basic measure of synchrony of spiking from all neurons.
:param spiketimes_set: Dictionary returned using :meth:`~analyseur.cbgtc.loader.LoadSpikeTimes.get_spiketimes_superset`
or using :meth:`~analyseur.cbgtc.loader.LoadSpikeTimes.get_spiketimes_subset`
:param binsz: integer or float; `0.01` [default]
:param window: Tuple in the form `(start_time, end_time)`; `(0, 10)` [default]
:param windowsz: integer or float, sliding window size; `0.5` [default]
:return: a number
NOTE: The computation is done on a sliding window resulting in a smoother frequency estimation otherwise
this is the same as :py:meth:`.compute_basic`. Therefore, unlike the simple histogram-based approach the use
of overlapping windows (sliding windows) results in a smoother frequency estimation.
.. raw:: html
<hr style="border: 2px solid red; margin: 20px 0;">
"""
# ============== DEFAULT Parameters ==============
if window is None:
window = cls.__siganal.window
if binsz is None:
binsz = cls.__siganal.binsz_100perbin
if windowsz is None:
windowsz = 0.5
#[spike_arrays, window] = cls.__get_spikearray_and_window(spiketimes_set, window, neurons="all")
spike_arrays = [np.array(neuron_spikes) for neuron_spikes in spiketimes_set.values()]
n_neurons = len(spike_arrays)
# Create evaluation time points
eval_times = np.arange(window[0] + windowsz/2, window[1] - windowsz, binsz)
n_times = len(eval_times)
freq_matrix = np.zeros((n_neurons, n_times))
for i, spike_times in enumerate(spike_arrays):
for t, t_center in enumerate(eval_times):
start = t_center - windowsz/2
stop = t_center + windowsz/2
spikes_in_slidingwindow = np.sum((spike_times >= start) & (spike_times <= stop))
freq_matrix[i, t] = spikes_in_slidingwindow / windowsz
# Remove time points with no activity
valid_points = np.sum(freq_matrix, axis=0) > 0
freq_matrix = freq_matrix[:, valid_points]
if freq_matrix.size == 0:
return 0.0, np.array([0]), np.array([0])
time_bins = np.arange(window[0], window[1] + binsz, binsz)
time_bins_center = (time_bins[:-1] + time_bins[1:]) / 2
[s_sync, _, _] = cls.__compute_for_basic(freq_matrix, filter_zeros=True)
return s_sync, freq_matrix, time_bins_center
[docs]
@classmethod
def compute_fano_factor(cls, spiketimes_set, binsz=None, window=None, filter_zeros=True):
"""
Returns the Fano factor as a measure of synchrony of spiking from all neurons.
:param spiketimes_set: Dictionary returned using :meth:`~analyseur.cbgtc.loader.LoadSpikeTimes.get_spiketimes_superset`
or using :meth:`~analyseur.cbgtc.loader.LoadSpikeTimes.get_spiketimes_subset`
:param binsz: integer or float; `0.01` [default]
:param window: Tuple in the form `(start_time, end_time)`; `(0, 10)` [default]
:param filter_zeros: boolean; `True` [default]
:return:
* fanofactor : float = variance / mean of population spike counts
* spike_matrix : numpy array, Population spike counts per time bin
* time_bins_center : numpy array, Center times of each bin
**Formula**
.. table:: Formula
====================================================================================================================== =====================================================================
Definitions Interpretation
====================================================================================================================== =====================================================================
total neurons, :math:`n_{Nuc}` total number of neurons in the Nucleus
neuron index, :math:`i` i-th neuron in the pool of :math:`n_{Nuc}` neurons
spike count, :math:`p^{(i)}(t)` spike count of the i-th neuron at time :math:`t`
count matrix, :math:`P = \\left[p(a,b) = p^{(a)}(b)\\right]_{\\forall{a \\in [1, n_{Nuc}], b \\in [t_0, t_T]}}` spike counts of all (:math:`n_{Nuc}`) neurons for all times
====================================================================================================================== =====================================================================
Let the :math:`var(\\cdot)`, `variance function <https://numpy.org/doc/stable/reference/generated/numpy.var.html>`_ and
the :math:`\\mu(\\cdot)`, `arithmetic mean function <https://numpy.org/doc/stable/reference/generated/numpy.mean.html>`_
be implemented as shown
.. math::
P = \\overset{\\begin{matrix}t_0 & \\quad\\quad & t_1 & & & &\\ldots & & & t_T\\end{matrix}}
{\\underset{
\\begin{matrix}
\\quad\\quad\\uparrow & \\quad\\quad\\quad & \\uparrow & \\quad &\\ldots & & & \\uparrow \n
\\quad\\Sigma_{t_0} & \\quad\\quad\\quad & \\Sigma_{t_1} & \\quad &\\ldots & & & \\Sigma_{t_T} & \\rightarrow var_{\\forall{t}} \n
\\quad\\Sigma_{t_0} & \\quad\\quad\\quad & \\Sigma_{t_1} & \\quad &\\ldots & & & \\Sigma_{t_T} & \\rightarrow \\mu_{\\forall{t}}
\\end{matrix}}
{\\begin{bmatrix}
p^{(1)}(t_0) & p^{(1)}(t_1) & \\ldots & p^{(1)}(t_T) \n
p^{(2)}(t_0) & p^{(2)}(t_1) & \\ldots & p^{(2)}(t_T) \n
\\vdots & \\vdots & \\ldots & \\vdots \n
p^{(i)}(t_0) & p^{(i)}(t_1) & \\ldots & p^{(i)}(t_T) \n
\\vdots & \\vdots & \\ldots & \\vdots \n
p^{(n_{Nuc})}(t_0) & p^{(n_{Nuc})}(t_1) & \\ldots & p^{(n_{Nuc})}(t_T)
\\end{bmatrix}
}}
We define
.. math::
A &\\triangleq var\\left(\\begin{bmatrix}
\\Sigma_{t_0} & \\Sigma_{t_1} & \\ldots & \\Sigma_{t_T}
\\end{bmatrix}\\right) = var_{\\forall{t}} \n
B &\\triangleq \\mu\\left(\\begin{bmatrix}
\\Sigma_{t_0} & \\Sigma_{t_1} & \\ldots & \\Sigma_{t_T}
\\end{bmatrix}\\right) = \\mu_{\\forall{t}}
Then, synchrony is measured as
.. math::
F_{Sync} = \\frac{A}{B} = \\frac{var\\left(\\left[\\sum_{\\forall{i}}p^{(i)}(t)\\right]_{\\forall{t}}\\right)}{\\mu\\left(\\left[\\sum_{\\forall{i}}p^{(i)}(t)\\right]_{\\forall{t}}\\right)}
**INTERPRETATION**
.. table:: Interpretation
================= ===================================================
Fano Factor Interpretation
================= ===================================================
:math:`= 1` Poisson-like variability
:math:`< 1` More regular/periodic than Poisson (Sub-Poisson)
:math:`> 1` More variable/bursty than Poisson (Super-Poisson)
================= ===================================================
COMMENTS:
* Fano factor can indicate population-level synchrony
- higher values correspond to more synchrony
* Fano factor depends on time bin size.
**Notes on Filtering**
The default argument `filter_zeros=True` means only spike counts during active periods are accounted for
making the computation. This is the default because,
* Biological relevance
- interest in variability of activity during firing periods **not during silent periods**.
* Statistical robustness
- zero bins artificially lower both variance and mean values.
* Practical interpretation
- Fano factor is more meaningful when computed over active periods.
.. raw:: html
<hr style="border: 2px solid red; margin: 20px 0;">
"""
#============== DEFAULT Parameters ==============
if window is None:
window = cls.__siganal.window
if binsz is None:
binsz = cls.__siganal.binsz_100perbin
# spike_matrix, time_bins = cls.__get_spike_matrix(spiketimes_set, window, binsz)
spike_matrix, _, time_bins = Rate.get_count_rate_matrix(spiketimes_set=spiketimes_set,
window=window, binsz=binsz,
neurons="all")
time_bins_center = (time_bins[:-1] + time_bins[1:]) / 2
# print(f"Unfiltered spike matrix: {spike_matrix}")
[fanofactor, _] = cls.__compute_fano(spike_matrix, filter_zeros=filter_zeros)
return fanofactor, spike_matrix, time_bins_center
[docs]
@classmethod
def compute_fano_factor_multibins(cls, spiketimes_set, binsz_list=None, window=None):
"""
Returns the Fano factor for multiple bin sizes
:param spiketimes_set: Dictionary returned using :meth:`~analyseur.cbgtc.loader.LoadSpikeTimes.get_spiketimes_superset`
or using :meth:`~analyseur.cbgtc.loader.LoadSpikeTimes.get_spiketimes_subset`
:param binsz_list: list of integer or float
:param window: Tuple in the form `(start_time, end_time)`; `(0, 10)` [default]
:return: For each binsz in binsz_list
* fanofactor : float = variance / mean of population spike counts
* spike_matrix : numpy array, Population spike counts per time bin
* time_bins_center : numpy array, Center times of each bin
**Formula**
see :py:meth:`.compute_fano_factor`
.. raw:: html
<hr style="border: 2px solid red; margin: 20px 0;">
"""
# ============== DEFAULT Parameters ==============
if window is None:
window = cls.__siganal.window
if binsz_list is None:
binsz_list = [
cls.__siganal.binsz_10perbin,
cls.__siganal.binsz_100perbin,
cls.__siganal.binsz_1000perbin
]
results = {}
for binsz in binsz_list:
# spike_matrix, time_bins = cls.__get_spike_matrix(spiketimes_set, window, binsz)
spike_matrix, _, time_bins = Rate.get_count_rate_matrix(spiketimes_set=spiketimes_set,
window=window, binsz=binsz,
neurons="all")
time_bins_center = (time_bins[:-1] + time_bins[1:]) / 2
[fanofactor, _] = cls.__compute_fano(spike_matrix)
results[binsz] = {
"fanofactor": fanofactor,
"spike_matrix": spike_matrix,
"time_bins_center": time_bins_center,
# variance of population spike counts during active periods, i.e, excluding time bins with no activity
# "variance": np.var(spike_matrix[spike_matrix > 0]) if np.any(spike_matrix > 0) else 0,
# mean of population spike counts during active periods, i.e, excluding time bins with no activity
# "mean": np.mean(spike_matrix[spike_matrix > 0]) if np.any(spike_matrix > 0) else 0,
}
return results
[docs]
@classmethod
def compute_sync_index_from_PSTH(cls, spiketimes_set, binsz=None, window=None):
"""
Synchrony as variance of total population rate
.. raw:: html
<hr style="border: 2px solid red; margin: 20px 0;">
"""
[counts, bin_info, _, _, _] = PSTH.compute_poolPSTH(spiketimes_set, neurons="all",
binsz=binsz, window=window)
# Total Population Rate
pop_rate = counts / bin_info["binsz"]
# Synchrony as variance of total population rate
sync_index = np.var(pop_rate)
return sync_index
[docs]
@classmethod
def compute_pairwise_ci(cls, spiketimes_set, binsz=None, window=None):
"""
Calculate average pairwise correlation index.
.. raw:: html
<hr style="border: 2px solid red; margin: 20px 0;">
"""
# ============== DEFAULT Parameters ==============
if window is None:
window = cls.__siganal.window
if binsz is None:
binsz = cls.__siganal.binsz_100perbin
[desired_spiketimes_subset, _] = get_desired_spiketimes_subset(spiketimes_set, neurons="all")
n_neurons = len(desired_spiketimes_subset)
duration = window[1] - window[0]
pairwise_ci = []
for i in range(n_neurons):
for j in range(i+1, n_neurons):
# Calculate firing rates
r_i = len(desired_spiketimes_subset[i]) / duration
r_j = len(desired_spiketimes_subset[j]) / duration
# Count coincidences
coincidences = 0
for spike_i in desired_spiketimes_subset[i]:
for spike_j in desired_spiketimes_subset[j]:
if abs(spike_i - spike_j) <= binsz:
coincidences += 1
break # count each spike only once
# Expected coincidences by chance
expected = 2 * binsz * duration * r_i * r_j
# Correlation Index
if expected > 0:
ci = (coincidences - expected) / expected
else:
ci = 0
pairwise_ci.append(ci)
return np.mean(pairwise_ci) if pairwise_ci else 0
[docs]
@classmethod
def compute_pairwise_corr(cls, spiketimes_set, binsz=None, window=None):
"""
Returns the Pairwise correlation as a measure of synchrony of spiking from all neurons.
:param spiketimes_set: Dictionary returned using :meth:`~analyseur.cbgtc.loader.LoadSpikeTimes.get_spiketimes_superset`
or using :meth:`~analyseur.cbgtc.loader.LoadSpikeTimes.get_spiketimes_subset`
:param binsz: integer or float; `0.01` [default]
:param window: Tuple in the form `(start_time, end_time)`; `(0, 10)` [default]
:return: a number
**Formula**
.. table:: Formula
====================================================================================================================== =====================================================================
Definitions Interpretation
====================================================================================================================== =====================================================================
total neurons, :math:`n_{Nuc}` total number of neurons in the Nucleus
neuron index, :math:`i` i-th neuron in the pool of :math:`n_{Nuc}` neurons
spike count, :math:`p^{(i)}(t)` spike count of the i-th neuron at time :math:`t`
count matrix, :math:`P = \\left[p(a,b) = p^{(a)}(b)\\right]_{\\forall{a \\in [1, n_{Nuc}], b \\in [t_0, t_T]}}` spike counts of all (:math:`n_{Nuc}`) neurons for all times
====================================================================================================================== =====================================================================
Let the :math:`var(\\cdot)`, `variance function <https://numpy.org/doc/stable/reference/generated/numpy.var.html>`_ and
the :math:`\\mu(\\cdot)`, `arithmetic mean function <https://numpy.org/doc/stable/reference/generated/numpy.mean.html>`_
be implemented as shown
.. math::
P = \\overset{\\begin{matrix}t_0 & \\quad\\quad & t_1 & & & &\\ldots & & & t_T\\end{matrix}}
{\\underset{
\\begin{matrix}
\\quad\\quad\\uparrow & \\quad\\quad\\quad & \\uparrow & \\quad &\\ldots & & & \\uparrow \n
\\quad\\pi_{t_0} & \\quad\\quad\\quad & \\pi_{t_1} & \\quad &\\ldots & & & \\pi_{t_T} & \\rightarrow var_{\\forall{t}} \n
\\quad\\pi_{t_0} & \\quad\\quad\\quad & \\pi_{t_1} & \\quad &\\ldots & & & \\pi_{t_T} & \\rightarrow \\mu_{\\forall{t}}
\\end{matrix}}
{\\begin{bmatrix}
p^{(1)}(t_0) & p^{(1)}(t_1) & \\ldots & p^{(1)}(t_T) \n
p^{(2)}(t_0) & p^{(2)}(t_1) & \\ldots & p^{(2)}(t_T) \n
\\vdots & \\vdots & \\ldots & \\vdots \n
p^{(i)}(t_0) & p^{(i)}(t_1) & \\ldots & p^{(i)}(t_T) \n
\\vdots & \\vdots & \\ldots & \\vdots \n
p^{(n_{Nuc})}(t_0) & p^{(n_{Nuc})}(t_1) & \\ldots & p^{(n_{Nuc})}(t_T)
\\end{bmatrix}
}}
We define
.. math::
A &\\triangleq var\\left(\\begin{bmatrix}
\\pi_{t_0} & \\pi_{t_1} & \\ldots & \\pi_{t_T}
\\end{bmatrix}\\right) = var_{\\forall{t}} \n
B &\\triangleq \\mu\\left(\\begin{bmatrix}
\\pi_{t_0} & \\pi_{t_1} & \\ldots & \\pi_{t_T}
\\end{bmatrix}\\right) = \\mu_{\\forall{t}}
Then, synchrony is measured as
.. math::
F_{Sync} = \\frac{A}{B} = \\frac{var\\left(\\left[\\sum_{\\forall{i}}p^{(i)}(t)\\right]_{\\forall{t}}\\right)}{\\mu\\left(\\left[\\sum_{\\forall{i}}p^{(i)}(t)\\right]_{\\forall{t}}\\right)}
.. raw:: html
<hr style="border: 2px solid red; margin: 20px 0;">
"""
#============== DEFAULT Parameters ==============
if window is None:
window = cls.__siganal.window
if binsz is None:
binsz = cls.__siganal.binsz_100perbin
# spike_matrix, _, time_bins = Rate.get_count_rate_matrix(spiketimes_set=spiketimes_set,
# window=window, binsz=binsz,
# neurons="all")
spike_matrix, active_neurons, time_bins = cls.__count_matrix_and_active_neurons(spiketimes_set, window, binsz)
time_bins_center = (time_bins[:-1] + time_bins[1:]) / 2
[correlations, pairs] = cls.__compute_pairwise_corr(spike_matrix, active_neurons)
return correlations, pairs, spike_matrix, time_bins_center