Spiking data analysis & processing (spynal.spikes)

Preprocessing, basic analyses, and plotting of neural spiking activity

Overview

Functionality includes computing spike rates (using binning or spike density methods) and their statistict, inter-spike intervals and their statistics, and spike data preprocessing and plotting.

Most functions expect one of two formats of spiking data:

• bool

  Binary spike trains where 1’s label times of spikes and 0’s = no spike, in a Numpy ndarray of dtype=bool, where one axis corresponds to time, but otherwise can have any arbitrary dimensionality (including other optional axes corresponding to trials, units, etc.)

• timestamp

  Explicit spike timestamps in a Numpy ndarray of dtype=object. This is a “ragged nested sequence” (array of lists/sub-arrays of different length), analogous to Matlab cell arrays. Each object element is a variable-length list-like 1D subarray of spike timestamps (of dtype float or int), and the container object array can have any arbitrary dimensionality (including optional axes corresponding to trials, units, etc.).

  Alternatively, for a single spike train, a simple 1D list or ndarray of timestamps may be given instead.

Most functions perform operations in a mass-univariate manner. This means that rather than embedding function calls in for loops over units, trials, etc., like this:

for unit in units:
    for trial in trials:
        results[trial,unit] = compute_something(data[trial,unit])

You can instead execute a single call on ALL the data, labeling the relevant axis for the computation (usually time here), and it will run in parallel (vectorized) across all units, trials, etc. in the data, like this:

results = compute_something(data, axis)

Function list

Spike count/rate and inter-spike interval computation

• rate : Estimate spike rates/counts using given method

• bin_rate : Compute spike counts/rates in series of time bins (regular or not)

• density : Compute spike density (smoothed rate) with given kernel

• isi : Compute inter-spike intervals from spike data

Rate and inter-spike interval stats

• rate_stats : Compute given statistic on spike rate data

• isi_stats : Compute given statistic on inter-spike interval data

Spike waveform-shape stats

• waveform_stats : Compute given statistic on spike waveform data

• width : Trough-to-peak temporal width of waveform

• repolarization : Time for waveform to decay after peak

• trough_width : Width of waveform trough

• amp_ratio : Ratio of trough/peak amplitude

Preprocessing

• times_to_bool : Convert spike timestamps to binary spike trains

• bool_to_times : Convert binary spike train to timestamps

• cut_trials : Cut spiking data into trials

• select_time_range : Select specific time range in spike data, deleting data outside range

• realign_data : Realign data to new within-trial times (new t=0)

• pool_electrode_units : Pool all units on each electrode into a multi-unit

Plotting

• plot_raster : Generate a raster plot

• plot_mean_waveforms : Plot mean spike waveforms from one/more units

• plot_waveform_heatmap : Plot heatmap (2d histogram) of spike waveforms

Synthetic data generation

• simulate_spike_rates : Generate synthetic rates from Poisson distribution

• simulate_spike_trains : Generate synthetic spike trains from Poisson process

Function reference

rate(data, method='bin', **kwargs)

Estimate spike rates (or counts) using given method

Spiking data can be timestamps or binary (0/1) spike trains

Parameters:

• data (ndarray, shape=Any, dtype=object (each element = (n_spikes,) array) or) –

  ndarray, shape=(…,n_timepts,…), dtype=bool

  List(s) of spike timestamps (in s). Can be given for either a single spike train, or for multiple spike trains (eg different trials, units, etc.) within an object array of any arbitrary shape. -or- Binary/boolean representation of spike times, for either a single or multiple spike trains.

• method ({'bin','density'}, default: 'bin') –

  Spike rate estimation method:

  • ’bin’ : Traditional rectangular-binned rate (aka PSTH; see bin_rate())

  • ’density’ : Kernel density estimator for smoothed spike rate (see density())

• **kwargs – Any further arguments passed as-is to specific rate estimation function.

Returns:

• rates (ndarray, shape=(…,n_timepts_out) or (…,n_timepts_out,…)) – Estimated spike rates in spk/s (or spike counts) using given method (and for each trial/unit/etc. in data).

  For timestamp data, same dimensionality as data, with time axis appended to end. For boolean data, same dimensionality as data. If only a single spike train is input, output is (n_timepts_out,) vector.

  For ‘bin’ method, n_timepts_out = n_bins. For ‘density’ method, n_timepts_out depends on lims, smp_rate

  dtype is float (except int for method == ‘bin AND output = ‘count’)

• timepts (ndarray, shape=(n_timepts_out,) or (n_bins,2)) – For ‘density’ method: Time sampling vector (in s). shape=(n_timepts_out,). For ‘bin’ method: [start,end] time of each time bin (in s). shape=(n_bins,2).

bin_rate(data, lims=None, width=0.05, step=None, bins=None, output='rate', axis=-1, timepts=None)

Compute spike rate/count within given sequence of hard-edged time bins

Spiking data can be timestamps or binary (0/1) spike trains

Use lims/width/step to set standard-width sliding window bins or use bins to set any arbitrary custom time bins

NOTE: Spikes are counted within each bin including the start, but excluding the end of the bin. That is each bin is defined as the half-open interval [start,end).

Parameters:

• data (ndarray, shape=Any, dtype=object (each element = (n_spikes,) array)) –

  or ndarray, shape=(…,n_timepts,…), dtype=bool

  List(s) of spike timestamps (in s). Can be given for either a single spike train, or for multiple spike trains (eg different trials, units, etc.) within an object array of any arbitrary shape. -or- Binary/boolean representation of spike times, for either a single or multiple spike trains.

• lims (array-like, shape=(2,)) – Full time range of analysis ([start,end] in s). Must input a value unless explicitly setting custom bins.

• width (scalar, default: 0.050 (50 ms)) – Full width (in s) of each time bin

• step (scalar, default: width (each bin starts at end of previous bin)) – Spacing (in s) between successive time bins

• bins (array-like, shape=(n_bins,2), default: setup_sliding_windows(width,lims,step)) – Alternative method for setting bins; overrides width/spacing/lims. [start,end] (in s) of each custom time bin. Bins can have any arbitrary width and spacing. Default generates bins with width and spacing from lims[0] to lims[1]

• output (str, default: 'rate') –

  Which type pf spike measure to return:

  • ’rate’ : spike rate in each bin, in spk/s. Float valued.

  • ’count’ : spike count in each bin. Integer valued.

  • ’bool’ : binary presence/absence of any spikes in each bin. Boolean valued.

• axis (int, default: -1 (last axis of array)) – Axis of binary data corresponding to time dimension. Not used for spike timestamp data.

• timepts (ndarray, shape=(n_timepts,)) – Time sampling vector (in s) for binary data. Not used for spike timestamp data, but MUST be input for binary data.

Returns:

• rates (ndarray, shape=(…,n_bins) or (…,n_bins,…), dtype=float or int or bool) – Spike rates (in spk/s) or spike counts in each time bin (and for each trial/unit/etc. in data).

  For timestamp data, same shape as data, with time axis appended to end. If only a single spike train is input, output is (n_bins,) vector.

  For boolean data, same shape as data with time axis length reduced to n_bins.

  dtype is float for output = ‘rate’, int for ‘count’, bool for ‘bool’.

• bins (ndarray, shape=(n_bins,2)) – [start,end] time (in s) of each time bin

psth(data, lims=None, width=0.05, step=None, bins=None, output='rate', axis=-1, timepts=None)

Alias of bin_rate(). See there for details.

density(data, lims=None, width=None, step=0.001, kernel='gaussian', buffer=None, axis=-1, timepts=None, **kwargs)

Compute spike density function (smoothed rate) via convolution with given kernel

Spiking data can be timestamps or binary (0/1) spike trains

NOTE: Spike densities exhibit “edge effects” where rate is biased downward at temporal limits of data, due to convolution kernel extending past edge of data and summing in zeros.

The best way to avoid this is to request lims that are well within the temporal limits of the data. We also ameliorate this by extending the requested lims temporarily by a buffer, which is trimmed off to the requested lims before returning the computed density.

Also, for binary spike data, if lims +/- buffer extend beyond the temporal limits of the data, we symmetrically reflect the data around its limits (note this can’t be done for timestamp data as we don’t know the data’s temporal limits).

Parameters:

• data (ndarray, shape=Any, dtype=object (each element = (n_spikes,) array)) –

  or ndarray, shape=(…,n_timepts,…), dtype=bool

  List(s) of spike timestamps (in s). Can be given for either a single spike train, or for multiple spike trains (eg different trials, units, etc.) within an object array of any arbitrary shape. -or- Binary/boolean representation of spike times, for either a single or multiple spike trains.

• lims (array-like, shape=(2,)) –

  Desired time range of returned spike density ([start,end] in s). For boolean spike data, defaults to (timepts[0],timepts[-1]). For spike timestamp data, a value MUST be input.

  NOTE: For either data type, lims should be within (or equal to) the full extent of the data itself. We test this for binary spike data, but we don’t know the actual extent of timestamp data, so that is up to the user!

• width (scalar, default: (kernel-dependent)) –

  Width parameter (in s) for given convolution kernel. Interpretation is specific to value of kernel:

  • ’gaussian’ : width = standard deviation, default: 0.50 (50 ms)

  • ’hanning’ : width = half-width (~ 2.53x Gaussian SD), default: 0.125 (125 ms)

• step (scalar, default: 0.001 (1 ms)) – Spacing (in s) between successive time points in returned spike density. Values > 1 ms imply some downsampling from original spike train time sampling. NOTE: Must be an integer multiple of original data sampling (1 ms for timestamp data). NOTE: This argument replaces downsmp in previous versions.

• kernel (str or ndarray or callable, default: 'gaussian') –

  Convolution kernel to use. Can be given as kernel name, with current options:

  • ’gaussian’ : Gaussian kernel

  • ’hanning’ : Hanning kernel

  Alternatively, can input any arbitrary kernel as an array of float values or as a custom function that can take any additional kwargs as arguments and return an array.

• buffer (float, default: (kernel-dependent, approximates length of kernel's edge effects)) –

  Length (in s) of symmetric buffer to add to each end of time dimension (and trim off before returning) to avoid edge effects. Kernel defaults:

  • ’gaussian’ : 3*`width` (3x std dev)

  • ’hanning’ : width (half-width of kernel)

• axis (int, default: -1 (last axis of array)) – Axis of binary data corresponding to time dimension. Not used for spike timestamp data.

• timepts (ndarray, shape=(n_timepts,)) – Time sampling vector (in s) for binary data. Not used for spike timestamp data, but MUST be input for binary data.

• **kwargs – All other kwargs passed directly to kernel function

Returns:

• rates (ndarray, shape=(…,n_timepts_out) or (…,n_timepts_out,…), dtype=float) – Spike density function – smoothed spike rates (in spk/s) estimated at each timepoint (and for each trial/unit/etc. in data).

  For timestamp data, same dimensionality as data, with time axis appended to end. For boolean data, same dimensionality as data. If only a single spike train is input, output is (n_timepts_out,) vector. n_timepts_out = iarange(lims[0],lims[1],step)

• timepts_out (ndarray, shape=(n_timepts_out,)) – Time sampling vector (in s) for rates

isi(data, axis=-1, timepts=None)

Compute inter-spike intervals of one or more spike trains

Spiking data can be timestamps or binary (0/1) spike trains

Parameters:

• data (ndarray, shape=Any, dtype=object (each element = (n_spikes,) array)) – or ndarray, shape=(…,n_timepts,…), dtype=bool List(s) of spike timestamps (in s). Can be given for either a single spike train, or for multiple spike trains (eg different trials, units, etc.) within an object array of any arbitrary shape. -or- Binary/boolean representation of spike times, for either a single or multiple spike trains.

• axis (int, default: -1 (last axis of array)) – Axis of binary data corresponding to time dimension. Not used for spike timestamp data.

• timepts (ndarray, shape=(n_timepts,)) – Time sampling vector (in s) for binary data. Not used for spike timestamp data, but MUST be input for binary data.

Returns:

ISIs – Time intervals between each successive pair of spikes in data (in same time units as data). For boolean data, output is converted to timestamp-like config with time axis removed. For timestamp data, same shape as data (but with 1 fewer item per array cell). Array cells with <= 1 spike returned as empty arrays. If only a single spike train is input, output is (n_spikes-1,) vector.

Return type:

ndarray, shape=(n_spikes-1,) or ndarray, dtype=object (each elem = (n_spikes-1,) array)

interspike_interval(data, axis=-1, timepts=None)

Alias of isi(). See there for details.

rate_stats(rates, stat='Fano', axis=None, **kwargs)

Compute given statistic on spike rates of one or more spike trains

Input data must be spike rates, eg as computed using rate()

Stats may be along one/more array axes (eg trials) or across entire data array

Parameters:

• rates (ndarray, shape=(...,n_obs,...)) – Spike rate data. Arbitrary shape.

• stat ({'Fano','CV'}, default: 'Fano') –

  Rate statistic to compute. Options:

  • ’Fano’ : Fano factor = var(rate)/mean(rate), using fano()

  • ’CV’ : Coefficient of Variation = SD(rate)/mean(rate), using cv()

• axis (int, default: None) – Array axis to compute rate statistics along (usually corresponding to distict trials/observations). If None, computes statistic across entire array (analogous to np.mean/var).

• **kwargs – Any additional keyword args passed directly to statistic computation function

Returns:

stats – Rate statistic computed on data. For 1d data or axis=None, a single scalar value is returned. Otherwise, it’s an array w/ same shape as rates, but with axis reduced to length 1.

Return type:

float or ndarray, shape=(…,1,…)

isi_stats(ISIs, stat='Fano', axis='each', **kwargs)

Compute given statistic on inter-spike intervals of one or more spike trains

Input data must be inter-spike intervals, eg as computed using isi()

Can request data to be pooled along one/more axes (eg trials) before stats computation

Parameters:

• ISIs (ndarray, shape=(n_spikes-1,) or ndarray, dtype=object (each elem = (n_spikes-1,) array)) – List of inter-spike intervals (in s), eg as computed by isi(). Can be given for either a single spike train, or for multiple spike trains (eg different trials, units, etc.) within an object array of any arbitrary shape.

• stat ({'Fano','CV','CV2','LV','burst_fract}, default: 'Fano') –

  ISI statistic to compute. Options:

  • ’Fano’ : Fano factor = var(ISIs)/mean(ISIs), using fano()

  • ’CV’ : Coefficient of Variation = SD(ISIs)/mean(ISIs), using cv()

  • ’CV2’ : Local Coefficient of Variation (Holt 1996), using cv2()

  • ’LV’ : Local Variation (Shinomoto 2009), using lv()

  • ’burst_fract’ : Measure of burstiness (% spikes in bursts), using burst_fract()

  CV2 and LV and CV-like measures that reduce influence of changes in spike rate on the metric by only measuring local variation (between temporally adjacent ISIs). See their specific functions for details.

• axis (int or None or 'each', default: 'each') –

  Axis of ISI data to pool ISIs along before computing stat. eg, for data that is shape (n_trials,n_units), if you want to compute a stat value for each unit, pooled across all trials, you’d set axis = 0. If axis=None, ISIs are pooled across the entire data array. If axis=’each’, stats are computed separately for each spike train in the array.

  NOTE: For locality-sensitive stats (‘CV2’,’LV’), axis MUST = ‘each’. NOTE: default ‘each’ is the opposite of default for rate_stats()

Returns:

stats – Given ISI stat, computed on ISI data. Return as a single scalar if axis=None. Otherwise returns as array of same shape as ISIs, but with axis reduced to length 1.

Return type:

float or ndarray

References

• Holt et al. (1996) Journal of Neurophysiology https://doi.org/10.1152/jn.1996.75.5.1806

• Shinomoto et al. (2009) PLoS Computational Biology https://doi.org/10.1371/journal.pcbi.1000433

burst_fract(ISIs, crit=0.02)

Compute measure of burstiness of ISIs of a spike train.

Burst_fract = fraction of all ISIs that are within a spike burst (ISI < 20 ms by default)

Parameters:

• ISIs (array-like, shape=(n_ISIs,)) – List of inter-spike intervals for a single spike train

• crit (float, default: 20 ms) – Criterion ISI value (s) to discriminate burst vs non-burst spikes

Returns:

burst – Fraction of all spikes that are within spike bursts

Return type:

float

waveform_stats(spike_waves, stat='width', axis=0, **kwargs)

Compute given statistic on one or more spike waveforms of one or more spike trains

Parameters:

• spike_waves (ndarray, shape=(...,n_timepts,...), dtype=float or) –

  ndarray, shape=Any, dtype=object (elem’s are (…,n_timepts,…) arrays) Spike waveform data, given in one of two formats:

  1. a single ndarray of waveform data with one or more waveforms. Shape is arbitrary, but axis should correspond to time samples of waveform(s).

  2. an object ndarray where each element contains waveform data like format (1) for one unit, trial, etc. Time axis and time sampling must be the same for all elements, but other dimensions need not be (ie there can be different numbers of spikes aacross trials, units, etc.)

• stat (str, default: 'width') –

  Spike waveform statistic to compute. Options:

  • ’width’Temporal width (s) btwn largest depol trough and subsequent hyperpol peak,

    using trough_to_peak_width()

  • ’repolarization’Time from hyperpol peak to subsequent inflection or percent decrease,

    using repolarization_time()

  • ’trough_width’ : Full with at half-height of depol trough, using trough_width()

  • ’amp_ratio’ : Ratio of amplitude of depol trough / hyperpol peak, using amp_ratio()

• axis (int, default: 0) – Axis of spike_waves to compute stat along, corresponding to waveform timepoints. For object array data ((2) above), this should be the axis of each object array element that corresponds to time.

Returns:

stats – ndarray, shape=Any, dtype=object (elem’s are (n_timepts,n_spikes) arrays) Given spike waveform stat, computed on each waveform in spike_waves. For 1d data (single waveform), a single scalar value is returned. Otherwise, it’s an array w/ same shape as spike_waves, but with axis reduced to length 1.

Return type:

float or ndarray, shape=(…,1,…) or

trough_to_peak_width(spike_wave, smp_rate)

Compute time difference between largest spike waveform trough (depolarization) and subsequent peak (after-hyperpolarization) for a single spike waveform

Parameters:

• spike_wave (ndarray, shape=(n_timepts,)) – Spike waveform to compute stat on

• smp_rate (float) – Sampling rate of spike waveform (Hz)

Returns:

stat – Temporal width of waveform (in s), from trough to after-hyperpolarization peak

Return type:

float

trough_width(spike_wave, smp_rate)

Compute full width at half-height of spike waveform trough (depolarization phase)

Parameters:

• spike_wave (ndarray, shape=(n_timepts,)) – Spike waveform to compute stat on

• smp_rate (float) – Sampling rate of spike waveform (Hz)

Returns:

stat – Temporal width (FWHM) of waveform depolarization trough (in s)

Return type:

float

repolarization_time(spike_wave, smp_rate, criterion=0.75)

Compute time of repolarization of after-hyperpolarization peak of single spike waveform – time difference from 1st peak after global trough to when it has decayed

Parameters:

• spike_wave (ndarray, shape=(n_timepts,)) – Spike waveform to compute stat on

• smp_rate (float) – Sampling rate of spike waveform data (Hz)

• criterion (float or 'inflection', default: 0.75) –

  What criterion is used to determine time of repolaration (decay to baseline) of peak:

  • floatUse the point where amplitude has decayed to given proportion of its peak value

    eg, for criterion=0.75, use point where amplitude is <= 75% of peak.

  • ’inflection’ : Use the inflection point (change in sign of the 2nd derivative)

Returns:

stat – Repolarization time (in s) of spike waveform

Return type:

float

trough_peak_amp_ratio(spike_wave)

Compute ratio of amplitudes (height) of largest spike waveform trough (depolarization) and subsequent peak (after-hyperpolarization) for a single spike waveform

Parameters:

spike_wave (ndarray, shape=(n_timepts,)) – Spike waveform to compute stat on

Returns:

stat – Amplitude ratio of waveform trough/peak

Return type:

float

bool_to_times(spike_bool, timepts, axis=-1)

Convert boolean (binary) spike train representaton to spike timestamps

Inverse function of times_to_bool()

Parameters:

• spike_bool (ndarray, shape=(...,n_timepts,...), dtype=bool) – Binary spike trains, where 1 indicates >= 1 spike in time bin, 0 indicates no spikes

• timepts (ndarray, shape=(n_timepts,)) – Time sampling vector for data (center of each time bin used to compute binary train)

• axis (int, default: -1 (last axis)) – Axis of data corresponding to time dimension

Returns:

spike_times – or ndarray, shape=(n_spikes,) Spike timestamps (in same time units as timepts), for each spike train in input. Returns as vector-valued array of timestamps if input is single spike train, otherwise as object array of variable-length timestamp vectors.

Return type:

ndarray, dtype=object (each element = (n_spikes,) array)

times_to_bool(spike_times, lims=None, width=0.001, bins=None)

Convert spike timestamps to boolean (binary) spike train representaton

Inverse function of bool_to_times()

Times bins for computing binary spike trains may be set either implicitly via lims and width, or set explicitly using bins.

Parameters:

• spike_times (ndarray, shape=Any, dtype=object (each element = (n_spikes,) array)) – or array-like, shape=(n_spikes,) List(s) of spike timestamps (in s). Can be given for either a single spike train, or for multiple spike trains (eg different trials, units, etc.) within an object array of any arbitrary shape.

• lims (array-like, shape=(2,)) – Full time range of analysis ([start,end] in s). Must input a value (unless explicitly setting custom bins)

• width (float, default: 0.001 (1 ms)) – Width of bin used to discretize spike times (s). Usually 1 ms.

• bins (array-like, (n_bins,2), default: setup_sliding_windows(width,lims,width)) – Alternative method for setting time bins. Overrides any values set for lims, width. [start,end] time (in s) of each custom time bin. Bins can in theory have any arbitrary width and spacing, but really you would always want equal width, and width = spacing, so each bin starts at end of last bin

Returns:

• spike_bool (ndarray, shape=(…,n_bins), dtype=bool) – Binary spike trains, where 1 indicates >= 1 spike in time bin, 0 indicates no spikes Same shape as spike_times, with time axis appended to end. If only a single spike train is input, output is (n_bins,) vector.

• timepts (ndarray, shape(n_bins,)) – Time sampling vector (in s). Center of each time bin used to compute binary spike data

cut_trials(data, trial_lims, smp_rate=None, axis=None, trial_refs=None)

Cut time-continuous spiking data into trials

Spiking data may be in form of spike timestamps or binary spike trains

Parameters:

• data (ndarray, shape=Any, dtype=object (elems=(n_spikes,) arrays)) –

  or ndarray, shape=(…,n_timepts,…), dtype=bool

  Time-continuous spiking data (not cut into trials). Arbitrary shape, could include multiple channels, etc. Given in one of two formats:

  • timestamp: Spike timestamps, usually in seconds referenced to some within-trial event

  • bool: Binary (1/0) spike trains in bool array

  For binary spike data, additional smp_rate and axis keyword arguments must be input to indicate the sampling rate (in Hz) and the array time axis.

• trial_lims (array-like, shape=(n_trials,2)) – List of [start,end] times of each trial (in same timebase as data) to use to cut data

• smp_rate (scalar) – Sampling rate of binary spiking data (Hz). Must input a value for binary spiking data; not used for spike timestamp data.

• axis (int, default: 0) – Axis of data array corresponding to time samples. Only used for binary spike data.

• trial_refs (array-like, shape=(n_trials,), default: None) – List giving event time in each trial to re-reference trial’s spike timestamps to (ie this sets t=0 for each trial). If None, just leave timestamps in original timebase). Only used for spike timestamp data.

Returns:

cut_data – or ndarray, shape=(…,n_trial_timepts,…,n_trials), dtype=bool Spiking data segmented into trials. Trial axis is appended to end of all axes in input data. Shape is otherwise the same for timestamp data, but for binary data time axis is reduced to length implied by trial_lims.

Return type:

ndarray, shape=(…,n_trials), dtype=object (elems=(n_trial_spikes,) arrays

select_time_range(data, time_range, time_axis=None, timepts=None)

Extract spiking data from given time range, deleting any spikes/data outside this range

For example, you may have spiking data extending over several seconds per trial, but only want spikes from t=0 to 1 s for analysis

Spiking data may be in form of spike timestamps or binary spike trains

Parameters:

• data (ndarray, shape=Any, dtype=object (elems=(n_spikes[trial,unit,etc.],) arrays)) –

  or ndarray, shape=(…,n_timepts,…), dtype=bool

  Spiking data, given in one of two formats:

  • timestamp: Spike timestamps, usually in seconds referenced to some within-trial event

  • bool: Binary (1/0) spike trains in bool array

  Can be any arbitrary shape (including having multiple units), as long as time_axis is given for bool data

• time_range (array-like, shape=(2,)) – Time range to select spiking data from, inclusive of endpoints. ie, all data such that time_range[0] <= data or <= time_range[1] is retained in output, and all other data is deleted.

• time_axis (int) – Axis of bool data corresponding to time samples. Must input for bool data; not used for spike timestamps.

• timepts (array-like, shape(n_timepts)) – Time sampling vector for bool data (in s). Must input for bool data; not used for spike timestamps.

Returns:

• sel_data (ndarray, shape=Any, dtype=object (elems=(n_spikes[trial,unit,etc.],) arrays)) – or ndarray, shape=(…,n_timepts_out,…), dtype=bool Data with only spikes within given time range retained. For timestamp data, this has same shape as input data. For binary data, time axis is reduced to length implied by time_range, but otherwise array has same shape as input data.

• sel_bool (ndarray, shape=Any, dtype=object (elems=(n_spikes[trial,unit,etc.],),dtype=bool arrays)) – or ndarray, shape=(n_timepts,) dtype=bool Boolean vector(s) used to select spikes. Useful for selecting same spikes out of a parallel data structure, such as spike waveforms with the same trials/units/etc. as spike timestamp data. For timestamp data, this has same shape as input data (but each object array element has dtype bool). For binary data, this has same shape as input data.

realign_data(data, align_times, trial_axis, time_axis=None, timepts=None, time_range=None)

Realign timing of trial-cut spiking data on new within-trial event times (eg new trial event) so that t=0 on each trial at given event.

For example, data aligned to a start-of-trial event might need to be realigned to the behavioral response.

Spiking data may be in form of spike timestamps or binary spike trains

Parameters:

• data (ndarray, shape=(...,n_trials,...), dtype=object (elems=(n_spikes[trial],) arrays)) –

  or ndarray, shape=(…,n_timepts,…), dtype=bool

  Spiking data, given in one of two formats:

  • timestamp: Spike timestamps, usually in seconds referenced to some within-trial event

  • bool: Binary (1/0) spike trains in bool array

  Can be any arbitrary shape (including having multiple units), as long as trial_axis is given (and also time_axis for bool data)

• align_times (array-like, shape=(n_trials,)) – New set of times (in old reference frame) to realign spiking data to

• trial_axis (int) – Axis of data corresponding to trials. Must input for either data type.

• time_axis (int) – Axis of bool data corresponding to time samples. Must input for bool data; not used for spike timestamps.

• timepts (array-like, shape(n_timepts)) – Time sampling vector for bool data (in s). Must input for bool data; not used for spike timestamps.

• time_range (array-like, shape=(2,)) – Time range to extract from each trial of bool data around new align time ([start,end] time in s relative to align_times). eg, time_range=(-1,1) -> extract 1 s on either side of align event. Must input for bool data; not used for spike timestamps.

Returns:

realigned – or ndarray, shape=(…,n_timepts_out,…), dtype=bool Data realigned to given within-trial times. For timestamp data, this has same shape as input data. For binary data, time axis is reduced to length implied by time_range, but otherwise array has same shape as input data.

Return type:

ndarray, shape=(…,n_trials,…), dtype=object (elems=(n_spikes[trial],) arrays)

realign_data_on_event(data, event_data, event, **kwargs)

Convenience wrapper around realign_data for relaligning to a given named event within a per-trial dataframe or dict variable.

Only parameters differing from realign_data() are described here.

Parameters:

• event_data (dict, {str : ndarray, shape=(n_trials,)} or DataFrame, shape=(n_trials,n_events)) – Per-trial event timing data to use to realign spike timestamps.

• event (str) – Dict key or DataFrame column name whose associated values are to be used to realign data

pool_electrode_units(data_sua, electrodes, axis=-1, elec_set=None, return_idxs=False, sort=True)

Pool spiking data across all units (neurons) on each electrode into a single multi-unit for each electrode

Spiking data may be in form of spike timestamps, binary spike trains, or spike rates/counts

Parameters:

• data_sua (ndarray, shape=(...,n_units,...), dtype=object (elems= (n_spikes[unit],) arrays)) –

  or dtype=bool or dtype=float or int

  Spiking data for multiple single units on one or more electrodes, and optionally for different trials/conditions/etc. (arbitrary shape), in 1 of 3 formats:

  • timestamp: Lists of spike timestamps (dtype=object)

  • bool: Binary (0/1) spike trains (dtype=bool)

  • rate: Sets of spike rates or counts (dtype=float or int)

• electrodes (array-like, shape=(n_units,)) – List of electrode numbers of each unit in data_sua

• axis (int, default: -1 (last axis)) – Axis of data_sua corresponding to different units.

• elec_set (array-like, shape=(n_elecs,), default: unsorted_unique(electrodes)) – Set of unique electrodes in electrodes. Default uses all unique values in electrodes.

• return_idxs (bool, default: False (only return data_mua)) – If True, additionally returns list of indexes corresponding to 1st occurrence of each electrode in electrodes

• sort (bool, default: True) – If True, sort pooled timestamps so they remain in sequential order after concatenation. Only used for timestamp data.

Returns:

• data_mua (ndarray, shape=(…,n_elecs,…), dtype=object (elems= (n_spikes[elec],) arrays)) – or dtype=bool or dtype=float or int Spiking data pooled across all single units into a single electrode-level multi-unit for each electrode. Same shape as data_sua, but with axis reduced to length=n_elecs.

  For timestamp data, spike timestamps are concatenated together across units, and optionally resorted into sequential order (if sort is True).

  For bool data, spike trains are combined via “logical OR” across units (spike is registered in output at times when there is a spike in any single-unit on electrode).

  For rate data, spike rates/counts are summed across units.

• elec_idxs (ndarray, shape=(n_elecs,), dtype=int, optional) – Indexes of 1st occurrence of each electrode in elec_set within electrodes. Can be used to transform any corresponding metadata appropriately. Only returned if return_idxs is True.

Examples

data_mua = pool_electrode_units(data_sua, electrodes, return_idxs=False)

data_mua, elec_idxs = pool_electrode_units(data_sua, electrodes, return_idxs=True)

plot_raster(data, ax=None, graphics=None, color='0.25', height=1.0, events=None, **kwargs)

Generate raster plot of spike times

Parameters:

• data (ndarray, shape=(n_spike_trains,), dtype=object (each element = (n_spikes,) array) or) –

  ndarray, shape=(n_spike_trains,n_timepts), dtype=bool

  List(s) of spike timestamps (in s). Can be given for either a single spike train, or for multiple spike trains (eg different trials, units, etc.) within an object array. -or- Binary/boolean representation of spike times, for either a single or multiple spike trains.

  Data is ultimately transformed to timestamp format for vector graphics plots, and to binary format for bitmap/raster graphics plots. Therefore, if the “wrong” format is input for a requested graphics type, you must also input additional keyword args required by the data format transformation functions – bool_to_times() requires a value timepts; times_to_bool()) requires a value for either lims or bins (see functions for details).

  NOTE: Since raster plots are 2-dimensional (time x trials/units/etc.), data here has stricter shape requirements than other functions – it can have at most 1 dimension other than time. That is, timestamp object arrays must be 1D (length=n_spike_trains); binary spike data must be <= 2D (shape=n_spike_trains,n_timepts).

• ax (Pyplot Axis object, default: plt.gca()) – Axis to plot into

• graphics ({'bitmap','vector'}, default: 'bitmap' for binary data, 'vector' for timestamps) –

  Type of graphics to use for plotting rasters:

  • Vector graphics using plt.plot(). Individual spikes rendered as line objects

  that are manipulatable in vector graphics software like Adobe Illustrator, but much slower.

  • Bitmapped (aka raster) graphics using :func:plt.imshow. Full raster plot rendered as

  single (unmanipulatable) image, but plots much faster. (Note: “raster graphics” is a term from computer graphics that is unrelated to “raster plots” in electrophysiology.)

  NOTE: In some cases, a substantial number of spikes “disappear” from bitmap plots. This seems to happen when plotting into smaller figures (eg inline/within-IDE plots or default-size separate-window figures), and when these figures are save to bitmap-based file formats (eg png, jpg). This problem doesn’t seem to happen when you plot into full-screen separate-window figures (see plots.full_figure()), when figures are saved into vector-based file formats (eg pdf, eps), or when figures are saved to bitmap-based formats (eg png) using our version of savefig() plots.savefig().

• color (Matplotlib color specifier, default: '0.25' (dark gray)) – Color to plot all spikes in

• height (float, default: 1.0 (each spike height is full range for its row in raster)) – Height of each plotted spike (in fraction of distance btwn spike trains). Only used for vector graphics version (height always 1 for bitmap/raster graphics).

• events (callable or array-like, shape=(n_events,)) – List of event values (eg times) to plot as markers on x-axis -or- callable function that will just plot the event markers itself. See plots.plot_markers() for details.

Returns:

ax – Axis for plot

Return type:

Pyplot Axis object

References

https://en.wikipedia.org/wiki/Vector_graphics

plot_mean_waveforms(spike_waves, timepts=None, plot_sd=True, ax=None, **kwargs)

Plot mean spike waveform for each of one or more units

Parameters:

• spike_waves (ndarray, shape=(n_units,), dtype=object (elems=(n_timepts,n_spikes) arrays)) – Spike waveforms for one or more units

• timepts (array-like, shape=(n_timepts,), default: 0:n_timepts) – Common time sampling vector for each spike waveform

• plot_sd (bool, default: True) – If True, also plots standard deviation of waves as fill

• ax (Pyplot Axis object, default: plt.gca()) – Axis to plot into

• **kwargs –

  Any additional keyword args are interpreted as parameters of plt.axes() (settable Axes object attributes), plt.plot() (Line2D object attributes), or plt.fill() (Polygon object attributes), including the following (with given default values):

  xlimarray-like, shape=(2,), default: (timepts[0],timepts[-1])

  x-axis limits

  xticklabels,yticklabelsarray-like, default: [] (no labels)

  Labels for x/y ticks

Returns:

• lines (List of Line2D objects) – ax.plot output. Allows access to line properties of line.

• patches (List of Polygon objects) – ax.fill output. Allows access to patch properties of fill.

• ax (Axis object) – Axis plotted into.

plot_waveform_heatmap(spike_waves, timepts=None, ylim=None, n_ybins=20, ax=None, cmap='jet', **kwargs)

Plot heatmap (2D hist) of all spike waveforms across one or more units

Parameters:

• spike_waves (ndarray, shape=(n_units,), dtype=object (elem's are (n_timepts,n_spikes) arrays)) – Spike waveforms for one or more units

• timepts (array-like, shape=(n_timepts,), default: 0:n_timepts) – Common time sampling vector for each spike waveform

• ylim (array-like, shape=(2,), default: [min,max] of given waveforms) – [min,max] waveform amplitude for generating 2D histograms

• n_ybins (int, default: 20) – Number of histogram bins to use for y (amplitude) axis

• ax (Pyplot Axis object, default: plt.gca()) – Axis to plot into

• cmap (str or Colormap object, default: 'jet') – Colormap to plot heat map in

Returns:

ax – Axis plotted into

Return type:

Pyplot Axis object

simulate_spike_rates(gain=5.0, offset=5.0, n_conds=2, n_trials=1000, window=1.0, count=False, seed=None)

Simulate Poisson spike rates across multiple conditions/groups with given condition effect size

Parameters:

• gain (scalar or array-like, shape=(n_conds,), default: 5.0 (5 spk/s btwn-cond diff)) – Spike rate gain (in spk/s) for each condition, which sets effect size. If scalar, interpeted as spike rate difference between each successive conditions. If array, interpreted as specific spike rate gain over baseline for each condition. Set = 0 to simulate no expected difference between conditions.

• offset (scalar, default: 5.0 (5 spk/s baseline)) – Baseline rate added to condition effects

• n_conds (int, default: 2) – Number of distinct conditions/groups to simulate

• n_trials (int, default: 1000) – Number of trials/observations to simulate

• window (scalar, default: 1.0 s) – Time window (in s) to count simulated spikes over. Set = 1 if you want spike counts, rather than rates.

• count (bool, default: False (compute rates)) – If True, return integer-valued spike counts. If False, return float-valued rates.

• seed (int, default: None) – Random generator seed for repeatable results. Set=None for unseeded random numbers.

Returns:

• rates (ndarray, shape=(n_trials,), dtype=float or int) – Simulated Poisson spike rates

• labels (ndarray, shape=(n_trials,), dtype=int) – Condition/group labels for each trial. Sorted in group order to simplify visualization.

simulate_spike_trains(gain=5.0, offset=5.0, n_conds=2, n_trials=1000, time_range=1.0, refractory=0, seed=None, data_type='timestamp')

Simulate Poisson spike trains across multiple conditions/groups with given condition effect size

Parameters:

• gain (scalar or array-like, shape=(n_conds,), default: 5.0 (5 spk/s btwn-cond diff)) – Spike rate gain (in spk/s) for each condition, which sets effect size. If scalar, interpeted as spike rate difference between each successive conditions. If array, interpreted as specific spike rate gain over baseline for each condition. Set = 0 to simulate no expected difference between conditions.

• offset (scalar, default: 5.0 (5 spk/s baseline)) – Baseline rate added to condition effects

• n_conds (int, default: 2) – Number of distinct conditions/groups to simulate

• n_trials (int, default: 1000) – Number of trials/observations to simulate

• time_range (scalar, default: 1 s) – Full time range (in s) to simulate spike train over

• refractory (scalar, default: 0 (no refractory)) –

  Absolute refractory period in which a second spike cannot occur after another. Set=0 for proper Poisson process with no refractory.

  NOTE: currently implemented by simply deleting spikes < refractory period, which affects rates and is thus not optimal

• seed (int, default: None) – Random generator seed for repeatable results. Set=None for unseeded random numbers.

• data_type ({'timestamp','bool'}, default: 'timestamp') –

  Format of output spike trains:

  • ’timestamp’ : Spike timestamps in s relative to trial starts

  • ’bool’ : Binary (0/1) vectors flagging spike times

Returns:

• trains (ndarray, shape=(n_trials,), dtype=object or) – ndarray, shape=(n_trials,n_timepts), dtype=bool Simulated Poisson spike trains, returned either as list of timestamps relative to trial start or as binary vector for each trial (depending on data_type).

• labels (ndarray, shape=(n_trials,), dtype=int) – Condition/group labels for each trial. Sorted in group order to simplify visualization.

simulate_spike_waveforms(trough_time=0.0003, peak_time=0.00075, trough_amp=0.4, peak_amp=0.25, trough_width=0.0002, peak_width=0.0004, time_range=(-0.0002, 0.0014), smp_rate=30000.0, noise=0.01, n_spikes=100, seed=None)

Simulate set of spike waveforms with given shape and amplitude parameters

Parameters:

• trough_time (float, default: 0.5 ms, 1 ms) – Time of waveform depolarization trough, after-hyperpolarization peak (in s relative to t=0 at simulated trigger time)

• peak_time (float, default: 0.5 ms, 1 ms) – Time of waveform depolarization trough, after-hyperpolarization peak (in s relative to t=0 at simulated trigger time)

• trough_amp (float, default: 0.4, 0.25) – Absolute amplitude of waveform depolarization trough, after-hyperpolarization peak (in mV)

• peak_amp (float, default: 0.4, 0.25) – Absolute amplitude of waveform depolarization trough, after-hyperpolarization peak (in mV)

• trough_width (float, default: 0.25 ms, 0.5 ms) – Width (2*Gaussian SD in s) of waveform depolarization trough, after-hyperpolarization peak

• peak_width (float, default: 0.25 ms, 0.5 ms) – Width (2*Gaussian SD in s) of waveform depolarization trough, after-hyperpolarization peak

• time_range (array-like, shape=(2,), default: (-0.5 ms, 1.5 ms)) – Full time range that spike waveforms are sampled on ((start,end) relative to trigger in s)

• smp_rate (float, default: 30 kHz) – Sampling rate of spike waveforms (in Hz)

• noise (float, default: 0.01) – Amplitude of additive Gaussian noise (in mV)

• n_spikes (int, default: 100) – Total number of spike waveforms to simulate

• seed (int, default: None) – Random generator seed for repeatable results. Set=None for unseeded random numbers.

Returns:

• spike_waves (ndarray, shape=(n_timepts,n_spikes)) – Set of simulated spike waveforms, based on given parameters + noise

• timepts (ndarray, shape=(n_timepts)) – Time sampling vector for all simulated waveforms (in s)