Electrophysiology: trace analysis

The electrophysiology library also contains methods to analyze intracellular recordings. To import the electrophysiology library:

from brian.library.electrophysiology import *

There is a series of example scripts in the examples/electrophysiology folder. Currently, most methods are related to the analysis of spike shape.


You can low-pass filter a trace as follows:

v_lp=lowpass(v, tau)

where tau is the time constant (cut-off frequency 1/(2*pi*tau)) and v is the trace (a vector of values). By default, tau is in units of the timestep. Alternatively, one can specify the timestep:

v_lp=lowpass(v, tau, dt=0.1*ms)

Spike analysis

Detecting spikes

The following function returns the time indexes of spike peaks in a trace v:

peaks=spike_peaks(v, vc=-10*mV)

where vc is the voltage criterion (we consider that there is a spike when v>vc). The algorithm works as follows. First, we identify positive crossings of the voltage criterion. Then, after each positive crossing, we look for the first local maximum (that is, when the voltage first starts decreasing). The last spike is treated differently because the peak may occur after the end of the recording, in which case the last element is considered as the peak.

It is possible to omit the voltage criterion vc. In this case, it is guessed with the following (rather complex) function:


The idea of this algorithm is to look at the trace in phase space (v,dv/dt). In this space, spikes tend to circle around some area which contains no trajectory. It appears that, somewhere in the middle of these circles, there is a voltage vc for which trajectories are either increasing (dv>0, upstroke of a spike) or decreasing (dv<0, downstroke of a spike) but never still (dv=0). This means that a positive crossing of this voltage always leads to a spike. We identify this voltage by looking for the largest interval of voltages (v1,v2) for which there is no sign change of dv/dt (over two successive timesteps), and we set vc=(v1+v2)/2, the middle of this interval.

As this method is rather complex, it is strongly advised to manually check whether it gives reasonable results.

Voltage reset

The average voltage reset after a spike is calculated as the average first minimum after a spike, with the following function:

reset=reset_potential(v, peaks=None, full=False)

The time indexes of spike peaks can be given (this may save some computation time). With the full=True option, the standard deviation is also returned.

Spike threshold

There are 3 ways to measure the spike threshold. The first derivative method uses a threshold criterion on the first derivative dv/dt to identify spike onsets:

onsets=spike_onsets(v, criterion=None, vc=None)

where criterion is the derivative criterion and vc is the voltage criterion to detect spikes. Note that the criterion is in units of voltage per time step. First, the algorithm detects spike peaks. Then for each spike, we look for the last local maximum of dv/dt before the spike, which should be the inflexion point of the spike. Then we identify the last time before the inflexion point when dv/dt is smaller than the criterion. The function returns the time indexes of the onsets, not their values (which are v[onsets]). The derivative criterion may be automatically determined, using the following function:

criterion=find_onset_criterion(v, guess=0.1, vc=None)

where guess is an optional initial guess for the optimization method. The algorithm is simple: find the criterion that minimizes the variability of onsets.

There are two other methods to measure spike thresholds, but they do not always give very good results (perhaps the trace should be preliminary filtered):

onsets2=spike_onsets_dv2(v, vc=None)
onsets3=spike_onsets_dv3(v, vc=None)

The first one finds the maximum of the second derivative d2v/dt2, the second one finds the maximum of d3v/dt3. These are global maxima in each interspike interval (it could be that looking for the last local maximum gives better results).

The following function returns the depolarization slope preceding each spike as an array:

slopes=slope_threshold(v, onsets=None, T=None)

In this function, spike onset indexes are passed through the onset keyword. The depolarization slope is calculated by linear regression over the T time bins preceding each spike. The result is in units of the time bin.

In a similar way, the following function returns the average membrane potential preceding each spike as an array:

vm_threshold(v, onsets=None, T=None):

Spike shape

The following function returns the average spike duration, defined as the time from onset to reset (next voltage minimum):


The onsets can be passed to save computation time, with the onsets keyword. With the option full=True, the function returns: the mean time from onset to peak, the mean time from onset down to same value (note that this may not be meaningful for some neurons), mean time from onset to next minimum, and standard deviations for these 3 values.

Note: this function may change.

The following function returns the average spike-triggered voltage:

shape=spike_shape(v, onsets=None, before=100, after=100)

If onsets is unspecified, it is calculated with the spike_onsets function. Note that you can align spikes on other times, for example peaks. The arguments before and after specify the number of time steps before and after the triger times.

Note: this should not be specific to spikes, it’s a stimulus-triggered average.

Spike mask

It is often useful to discard spikes from the trace to analyse it. The following function returns an array of booleans which are True in spikes:

spike_mask(v, spikes=None, T=None)

The starting point of each spike (time bin) is given by the spikes variable (default: onsets), and T is the duration of each spike in time bins. This function can then be used to select the subthreshold trace or the spikes: