File: examples-electrophysiology_threshold_analysis.txt

package info (click to toggle)
brian 1.4.3-1
  • links: PTS, VCS
  • area: main
  • in suites: sid, stretch
  • size: 23,436 kB
  • sloc: python: 68,707; cpp: 29,040; ansic: 5,182; sh: 111; makefile: 61
file content (63 lines) | stat: -rw-r--r-- 1,802 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
.. currentmodule:: brian

.. index::
   pair: example usage; load
   pair: example usage; plot
   pair: example usage; subplot
   pair: example usage; show

.. _example-electrophysiology_threshold_analysis:

Example: threshold_analysis (electrophysiology)
===============================================

Analysis of spike threshold.

Loads a current clamp voltage trace, compensates (remove electrode voltage)
and analyses the spikes.

::

    from brian import *
    from brian.library.electrophysiology import *
    import numpy
    
    dt=.1*ms
    Vraw = numpy.load("trace.npy") # Raw current clamp trace
    I = numpy.load("current.npy")
    V, _ = Lp_compensate(I, Vraw, dt) # Electrode compensation
    
    # Peaks
    spike_criterion=find_spike_criterion(V)
    print "Spike detected when V exceeds",float(spike_criterion/mV),"mV"
    peaks=spike_peaks(V,vc=spike_criterion) # vc is optional
    
    # Onsets (= spike threshold)
    onsets=spike_onsets(V,criterion=3*dt,vc=spike_criterion) # Criterion: dV/dt>3 V/s
    
    # Spike-triggered average of V
    STA=spike_shape(V, onsets=onsets, before=100, after=100)
    
    print "Spike duration:",float(spike_duration(V,onsets=onsets)*dt/ms),"ms"
    print "Reset potential:",float(reset_potential(V,peaks=peaks)/mV),"mV"
    
    # Spike threshold statistics
    slope=slope_threshold(V,onsets=onsets,T=int(5*ms/dt))
    
    # Subthreshold trace
    subthreshold=-spike_mask(V)
    
    t=arange(len(V))*dt
    subplot(221)
    plot(t/ms,V/mV,'k')
    plot(t[peaks]/ms,V[peaks]/mV,".b")
    plot(t[onsets]/ms,V[onsets]/mV,".r")
    subplot(222)
    plot(((arange(len(STA))-100)*dt)/ms,STA/mV,'k')
    subplot(223)
    plot(t[subthreshold]/ms,V[subthreshold]/mV,'k')
    subplot(224)
    plot(slope/ms,V[onsets]/mV,'.')
    show()