Spike train analysis¶
ephysiopy.common.spikingcalcs.SpikeCalcsGeneric(spike_times: np.ndarray, cluster: int | None = None, event_ts: np.ndarray | None = None, event_window: np.ndarray = np.array([-0.5, 0.5]), pos_sample_rate: float = 50.0, sample_rate: float = 30000.0)
¶
Methods:
| Name | Description |
|---|---|
acorr |
Calculates the autocorrelogram of a spike train. |
contamination_percent |
Returns the contamination percentage of a spike train. |
get_ifr |
Returns the instantaneous firing rate of the cluster |
get_ifr_power_spectrum |
Returns the power spectrum of the instantaneous firing rate of a cell |
get_shuffled_ifr_sp_corr |
Returns an nShuffles x nSamples sized array of shuffled |
ifr_sp_corr |
Calculates the correlation between the instantaneous firing rate and |
mean_isi_range |
Calculates the mean of the autocorrelation from 0 to n seconds. |
psch |
Calculate the peri-stimulus count histogram of a cell's spiking |
psth |
Calculate the PSTH of event_ts against the spiking of a cell |
responds_to_stimulus |
Checks whether a cluster responds to a laser stimulus. |
shuffle_isis |
Shuffle the spike train ISIs and return a new spike train with |
smooth_spike_train |
Returns a spike train the same length as num pos samples that has been |
theta_band_max_freq |
Calculates the frequency with the maximum power in the theta band (6-12Hz) |
theta_mod_idx |
Calculates a theta modulation index of a spike train based on the cells |
theta_mod_idxV2 |
This is a simpler alternative to the theta_mod_idx method in that it |
theta_mod_idxV3 |
Another theta modulation index score this time based on the method used |
Attributes:
| Name | Type | Description |
|---|---|---|
n_spikes |
Returns the number of spikes in the cluster |
n_spikes
property
¶
Returns the number of spikes in the cluster
Returns:
| Type | Description |
|---|---|
int
|
The number of spikes in the cluster |
acorr(Trange: np.ndarray = np.array([-0.5, 0.5]), **kwargs) -> BinnedData
¶
Calculates the autocorrelogram of a spike train.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
Trange
|
ndarray
|
The range of times to calculate the autocorrelogram over (default is [-0.5, 0.5]). |
array([-0.5, 0.5])
|
**kwargs
|
dict
|
Additional keyword arguments. |
{}
|
Returns:
| Type | Description |
|---|---|
BinnedData
|
Container for the binned data. |
contamination_percent(**kwargs) -> tuple
¶
Returns the contamination percentage of a spike train.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
**kwargs
|
Passed into the contamination_percent function. |
{}
|
Returns:
| Type | Description |
|---|---|
tuple of float
|
Q - A measure of refractoriness. R - A second measure of refractoriness (kicks in for very low firing rates). |
get_ifr(spike_times: np.array, **kwargs) -> np.ndarray
¶
Returns the instantaneous firing rate of the cluster
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
spike_times
|
ndarray
|
The times in seconds at which the cluster fired. |
required |
n_samples
|
int
|
The number of samples to use in the calculation. Practically this should be the number of position samples in the recording. |
required |
Returns:
| Type | Description |
|---|---|
ndarray
|
The instantaneous firing rate of the cluster |
get_ifr_power_spectrum() -> tuple[np.ndarray, ...]
¶
Returns the power spectrum of the instantaneous firing rate of a cell
Used to calculate the theta_mod_idxV3 score above
Returns:
| Type | Description |
|---|---|
tuple of np.ndarray
|
The frequency and power of the instantaneous firing rate |
get_shuffled_ifr_sp_corr(ts: np.array, speed: np.array, nShuffles: int = 100, **kwargs)
¶
Returns an nShuffles x nSamples sized array of shuffled instantaneous firing rate x speed correlations
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
ts
|
ndarray
|
the times in seconds at which the cluster fired |
required |
speed
|
ndarray
|
the speed vector |
required |
nShuffles
|
int
|
the number of times to shuffle the timestamp vector 'ts' |
100
|
**kwargs
|
Passed into ifr_sp_corr |
{}
|
Returns:
| Type | Description |
|---|---|
ndarray
|
A nShuffles x nSamples sized array of the shuffled firing rate vs speed correlations. |
ifr_sp_corr(ts, speed, minSpeed=2.0, maxSpeed=40.0, sigma=3, nShuffles=100, **kwargs)
¶
Calculates the correlation between the instantaneous firing rate and speed.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
ts
|
The times in seconds at which the cluster fired. |
required | |
speed
|
ndarray
|
Instantaneous speed (nSamples lenght vector). |
required |
minSpeed
|
float
|
Speeds below this value are ignored. |
2.0
|
maxSpeed
|
float
|
Speeds above this value are ignored. |
40.0
|
sigma
|
int
|
The standard deviation of the gaussian used to smooth the spike train. |
3
|
nShuffles
|
int
|
The number of resamples to feed into the permutation test. |
100
|
**kwargs
|
method: how the significance of the speed vs firing rate correlation is calculated |
{}
|
Examples:
An example of how I was calculating this is:
rng = np.random.default_rng() method = stats.PermutationMethod(n_resamples=nShuffles, random_state=rng)
See Also
See scipy.stats.PermutationMethod.
mean_isi_range(isi_range: float) -> float
¶
Calculates the mean of the autocorrelation from 0 to n seconds. Used to help classify a neuron's type (principal, interneuron, etc).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
isi_range
|
int
|
The range in seconds to calculate the mean over. |
required |
Returns:
| Type | Description |
|---|---|
float
|
The mean of the autocorrelogram between 0 and n seconds. |
psch(bin_width_secs: float) -> np.ndarray
¶
Calculate the peri-stimulus count histogram of a cell's spiking against event times.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
bin_width_secs
|
float
|
The width of each bin in seconds. |
required |
Returns:
| Name | Type | Description |
|---|---|---|
result |
ndarray
|
Rows are counts of spikes per bin_width_secs. Size of columns ranges from self.event_window[0] to self.event_window[1] with bin_width_secs steps; so x is count, y is "event". |
psth() -> tuple[list, ...]
¶
Calculate the PSTH of event_ts against the spiking of a cell
Returns:
| Type | Description |
|---|---|
x, y : list
|
|
The list of time differences between the spikes of the cluster
|
|
and the events (x) and the trials (y)
|
|
responds_to_stimulus(threshold: float, min_contiguous: int, return_activity: bool = False, return_magnitude: bool = False, **kwargs) -> namedtuple
¶
Checks whether a cluster responds to a laser stimulus.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
threshold
|
float
|
The amount of activity the cluster needs to go beyond to be classified as a responder (1.5 = 50% more or less than the baseline activity). |
required |
min_contiguous
|
int
|
The number of contiguous samples in the post-stimulus period for which the cluster needs to be active beyond the threshold value to be classed as a responder. |
required |
return_activity
|
bool
|
Whether to return the mean reponse curve. |
False
|
return_magnitude
|
int
|
Whether to return the magnitude of the response. NB this is either +1 for excited or -1 for inhibited. |
False
|
Returns:
| Type | Description |
|---|---|
namedtuple
|
With named fields "responds" (bool), "normed_response_curve" (np.ndarray), "response_magnitude" (np.ndarray) |
shuffle_isis() -> np.ndarray
¶
Shuffle the spike train ISIs and return a new spike train with the same number of spikes but a different ISI distributedon. Useful for creating null distributions of spike metrics.
smooth_spike_train(npos, sigma=3.0, shuffle=None)
¶
Returns a spike train the same length as num pos samples that has been smoothed in time with a gaussian kernel M in width and standard deviation equal to sigma.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
npos
|
int
|
The number of position samples captured. |
required |
sigma
|
float
|
The standard deviation of the gaussian used to smooth the spike train. |
3.0
|
shuffle
|
int
|
The number of seconds to shift the spike train by. Default is None. |
None
|
Returns:
| Type | Description |
|---|---|
ndarray
|
The smoothed spike train. |
theta_band_max_freq(**kws)
¶
Calculates the frequency with the maximum power in the theta band (6-12Hz) of a spike train's autocorrelogram.
This function is used to look for differences in theta frequency in different running directions as per Blair. See Welday paper - https://doi.org/10.1523/jneurosci.0712-11.2011
Returns:
| Type | Description |
|---|---|
float
|
The frequency with the maximum power in the theta band. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If the input spike train is not valid. |
theta_mod_idx(**kwargs) -> float
¶
Calculates a theta modulation index of a spike train based on the cells autocorrelogram.
The difference of the mean power in the theta band (6-11 Hz) and the mean power in the 1-50 Hz band is divided by their sum to give a metric that lives between 0 and 1
Returns:
| Type | Description |
|---|---|
float
|
The difference of the values at the first peak and trough of the autocorrelogram. |
Notes
This is a fairly skewed metric with a distribution strongly biased to -1 (although more evenly distributed than theta_mod_idxV2 below)
theta_mod_idxV2() -> float
¶
This is a simpler alternative to the theta_mod_idx method in that it calculates the difference between the normalized temporal autocorrelogram at the trough between 50-70ms and the peak between 100-140ms over their sum (data is binned into 5ms bins)
Returns:
| Type | Description |
|---|---|
float
|
The difference of the values at the first peak and trough of the autocorrelogram. |
Notes
Measure used in Cacucci et al., 2004 and Kropff et al 2015
theta_mod_idxV3(**kwargs) -> float
¶
Another theta modulation index score this time based on the method used by Kornienko et al., (2024) (Kevin Allens lab) see https://doi.org/10.7554/eLife.35949.001
Uses the binned spike train instead of the autocorrelogram as the input to the periodogram function (they use pwelch in R; periodogram is a simplified call to welch in scipy.signal)
The resulting metric is similar to that in theta_mod_idx above except that the frequency bands compared to the theta band are narrower and exclusive of the theta band
Produces a fairly normally distributed score with a mean and median pretty close to 0
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
**kwargs
|
Passed into get_ifr_power_spectrum |
{}
|
Returns:
| Type | Description |
|---|---|
float
|
The difference of the values at the first peak and trough of the autocorrelogram. |
ephysiopy.common.spikingcalcs.get_burstiness(isi_matrix: np.ndarray, whiten: bool = False, plot_pcs: bool = False) -> np.ndarray
¶
Returns the burstiness of a waveform.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
isi_matrix
|
ndarray
|
A matrix of normalized interspike intervals (ISIs) for the neurons. Rows are neurons, columns are ISI time bins. |
required |
Returns:
| Type | Description |
|---|---|
ndarray
|
|
Notes
Algorithm:
1) The interspike intervals between 0 and 60ms were binned into 2ms bins, and the area of the histogram was normalised to 1 to produce a probability distribution histogram for each neuron
2) A principal components analysis (PCA) is performed on the matrix of the ISI probability distributions of all neurons
3) Neurons were then assigned to two clusters using a k-means clustering algorithm on the first three principal components
4) a linear discriminant analysis performed in MATLAB (‘classify’) was undertaken to determine the optimal linear discriminant (Fishers Linear Discriminant) i.e., the plane which best separated the two clusters in a three-dimensional scatter plot of the principal components.
Training on 80% of the data and testing on the remaining 20% resulted in a good separation of the two clusters.
5) A burstiness score was assigned to each neuron which was calculated by computing the shortest distance between the plotted point for each neuron in the three-dimensional cluster space (principal components 1,2 and 3), and the plane separating the two clusters (i.e., the optimal linear discriminant).
6) To ensure the distribution of these burstiness scores was bimodal, reflecting the presence of two classes of neuron (‘bursty’ versus ‘non-bursty’), probability density functions for Gaussian mixture models with between one and four underlying Gaussian curves were fitted and the fit of each compared using the Akaike information criterion (AIC)
7) Optionally plot the principal components and the centres of the kmeans results
ephysiopy.common.spikingcalcs.mahal(u, v)
¶
Returns the L-ratio and Isolation Distance measures calculated on the principal components of the energy in a spike matrix.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
u
|
ndarray
|
The first set of waveforms. |
required |
v
|
ndarray
|
The second set of waveforms. |
required |
Returns:
| Type | Description |
|---|---|
ndarray
|
The Mahalanobis distances. |
Raises:
| Type | Description |
|---|---|
Warning
|
If input size mismatch, too few rows, or complex inputs are detected. |
ephysiopy.common.spikingcalcs.cluster_quality(waveforms: np.ndarray = None, spike_clusters: np.ndarray = None, cluster_id: int = None, fet: int = 1)
¶
Returns the L-ratio and Isolation Distance measures calculated on the principal components of the energy in a spike matrix.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
waveforms
|
ndarray
|
The waveforms to be processed. If None, the function will return None. |
None
|
spike_clusters
|
ndarray
|
The spike clusters to be processed. |
None
|
cluster_id
|
int
|
The ID of the cluster to be processed. |
None
|
fet
|
int
|
The feature to be used in the PCA calculation (default is 1). |
1
|
Returns:
| Type | Description |
|---|---|
tuple
|
A tuple containing the L-ratio and Isolation Distance of the cluster. |
Raises:
| Type | Description |
|---|---|
Exception
|
If an error occurs during the calculation of the L-ratio or Isolation Distance. |
ephysiopy.common.spikingcalcs.xcorr(x1: np.ndarray, x2: np.ndarray | None = None, Trange: np.ndarray | list = np.array([-0.5, 0.5]), binsize: float = 0.001, normed=False, **kwargs) -> BinnedData
¶
Calculates the ISIs in x1 or x1 vs x2 within a given range.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
x1
|
ndarray
|
The times of the spikes emitted by the first cluster in seconds. |
required |
x2
|
ndarray
|
The times of the spikes emitted by the second cluster in seconds. If None, x1 is used. |
None
|
Trange
|
ndarray or list
|
Range of times to bin up in seconds (default is [-0.5, 0.5]). |
array([-0.5, 0.5])
|
binsize
|
float
|
The size of the bins in seconds (default is 0.001). |
0.001
|
normed
|
bool
|
Whether to divide the counts by the total number of spikes to give a probability (default is False). |
False
|
**kwargs
|
dict
|
Additional keyword arguments. |
{}
|
Returns:
| Type | Description |
|---|---|
BinnedData
|
A BinnedData object containing the binned data and the bin edges. |
ephysiopy.common.spikingcalcs.fit_smoothed_curve_to_xcorr(xc: BinnedData, **kwargs) -> BinnedData
¶
Idea is to smooth out the result of an auto- or cross-correlogram with a view to correlating the result with another auto- or cross-correlogram to see how similar two of these things are.
Check Brandon et al., 2011?2012?
ephysiopy.common.spikingcalcs.contamination_percent(x1: np.ndarray, x2: np.ndarray | None = None, **kwargs) -> tuple
¶
Computes the cross-correlogram between two sets of spikes and estimates how refractory the cross-correlogram is.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
x1
|
ndarray
|
The first set of spikes. |
required |
x2
|
ndarray
|
The second set of spikes. If None, x1 is used. |
None
|
**kwargs
|
dict
|
Additional keyword arguments that can be fed into xcorr. |
{}
|
Returns:
| Type | Description |
|---|---|
tuple
|
A tuple containing: - Q (float): A measure of refractoriness. - R (float): A second measure of refractoriness (kicks in for very low firing rates). |
Notes
Taken from KiloSorts ccg.m
The contamination metrics are calculated based on an analysis of the 'shoulders' of the cross-correlogram. Specifically, the spike counts in the ranges ±5-25ms and