Analysis and plotting

Most plotting should be done with the PyLab commands, all of which are loaded when you import Brian. See:

for help on PyLab. The scientific library Scipy is also automatically imported by the instruction from brian import *.

The most useful plotting instruction is the Pylab function plot. A typical use with Brian is:

plot(t/ms,vm/mV)

where t is a vector of times with units ms and vm is a vector of voltage values with units mV. To display the figures on the screen, the function show() must be called once (this should be the last line of your script), except when using IPython with the Pylab mode (ipython -pylab).

Brian currently defines just two plotting functions of its own, raster_plot() and hist_plot(). In addition, the StateMonitor object has a plot() method.

Raster plots

Spike trains recorded by a SpikeMonitor can be displayed as raster plots:

S=SpikeMonitor(group)
...
raster_plot(S)

Usual options of the plot command can also be passed to raster_plot(). One may also pass several spike monitors as arguments.

State variable plots

State values recorded by a StateMonitor can also be plotted as follows:

M = StateMonitor(group, 'V', record=[0,1,2])
...
M.plot()

Realtime plotting

Both raster_plot() and StateMonitor.plot() have real-time versions which update as the simulation runs, for example:

G = NeuronGroup(...)
spikemon = SpikeMonitor(G)
statemon = StateMonitor(G, 'V', record=range(5))
ion()
subplot(211)
raster_plot(spikemon, refresh=10*ms, showlast=200*ms)
subplot(212)
statemon.plot(refresh=10*ms, showlast=200*ms)
run(1*second)
ioff()
show()

The ion() and ioff() command activate and deactivate Pylab’s interactive plotting mode. The refresh parameter specifies how often (in simulation time) to refresh the plot - smaller values will slow down the simulation. The showlast option only plots the most recent values.

With some IDEs, you may need to do something like the following at the beginning of your script to make interactive mode work:

import matplotlib
matplotlib.use('WXAgg')

This is because the default graphical backend can sometimes interact badly with the IDE. Other options to try are GTKAgg, QTAgg, TkAgg.

Statistics

Here are a few functions to analyse first and second order statistical properties of spike trains, defined as ordered lists of spike times:

  • Firing rate: firing_rate(spikes) where spikes is a spike train (list of spike times).
  • Coefficient of variation: CV(spikes).
  • Cross-correlogram: correlogram(T1,T2,width=20*ms,bin=1*ms,T=None) returns the cross-correlogram of spike trains T1 and T2 with lag in [-width,width] and given bin size. T is the total duration (optional) and should be greater than the duration of T1 and T2. The result the rate of coincidences in each bin, returned as an array.
  • Autocorrelogram: autocorrelogram(T0,width=20*ms,bin=1*ms,T=None) is the same as correlogram(T0,T0,width=20*ms,bin=1*ms,T=None).
  • Cross-correlation function: CCF(T1,T2,width=20*ms,bin=1*ms,T=None) returns the cross-correlation function of T1 and T2, which is the same as the cross-correlogram divided by the bin size (which makes the result independent of the bin size).
  • Autocorrelation function: ACF(T0,width=20*ms,bin=1*ms,T=None), same as CCF(T0,T0,width=20*ms,bin=1*ms,T=None).
  • Cross-covariance function: CCVF(T1,T2,width=20*ms,bin=1*ms,T=None) is the cross-correlation function of T1 and T2 minus for the cross-correlation of independent spike trains with the same rates (product of rates).
  • Auto-covariance function: ACVF(T0,width=20*ms,bin=1*ms,T=None) is the same as CCVF(T0,T0,width=20*ms,bin=1*ms,T=None).
  • Total correlation coefficient: total_correlation(T1,T2,width=20*ms,T=None) is the integral of the cross-covariance function divided by the rate of T1, typically (but not always) between 0 and 1.
  • Vector strength: vector_strength(spikes,period) returns the vector strength of the given spike train with respect to the period. If each spike time with phase phi is represented by a vector with angle phi, then the vector strength is the length of the average vector. It equals 1 for spikes with constant phase and 0 for homogeneous phase distributions.
  • Gamma precision factor: gamma_factor(source, target, delta) returns the gamma precision factor between source and target trains, with precision delta.

These functions return NaN (not a number) when a spike train is empty.