File: examples-frompapers_Muller_et_al_2011.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 (124 lines) | stat: -rw-r--r-- 4,082 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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
.. currentmodule:: brian

.. index::
   pair: example usage; std
   pair: example usage; xlim
   pair: example usage; title
   pair: example usage; show
   pair: example usage; errorbar
   pair: example usage; sqrt
   pair: example usage; xlabel
   pair: example usage; linspace
   pair: example usage; ylabel
   pair: example usage; ylim
   pair: example usage; mean

.. _example-frompapers_Muller_et_al_2011:

Example: Muller_et_al_2011 (frompapers)
=======================================

Interplay of STDP and input oscillations
----------------------------------------
Figure 4 from:
Muller L, Brette R and Gutkin B (2011) Spike-timing dependent plasticity and 
feed-forward input oscillations produce precise and invariant spike phase-locking. Front. 
Comput. Neurosci. 5:45. doi: 10.3389/fncom.2011.00045

Description:
In this simulation, a group of IF neurons is given a tonic DC input and a tonic AC input.
The DC input is mediated by current injection (neurons.I, line 62), and the AC input is 
mediated by Poisson processes whose rate parameters are oscillating in time. Each neuron in 
the group is given a different DC input, ensuring a unique initial phase. After two seconds 
of simulation (to integrate out any initial transients), the STDP rule is turned on 
(ExponentialSTDP, line 68), and the population of neurons converges to the theoretically 
predicted fixed point. As there is some noise in the phase due to the random inputs, the 
simulation is averaged over trials (50 in Figure 4, though 10 trials should be fine for 
testing).

The trials run in parallel on all available processors (10 trials take about
2 minutes on a modern PC).

::

    
    ### IMPORTS
    from brian import *
    import multiprocessing
    
    ### PARAMETERS
    N=5000
    M=10
    taum=33*ms                          
    tau_pre=20*ms
    tau_post=tau_pre
    Ee=0*mV
    vt=-54*mV
    vr=-70*mV
    El=-70*mV
    taue=5*ms
    f=20*Hz
    theta_period = 1/f
    Rm=200*Mohm
    a = linspace(51,65,num=M)       
    weights = .001
    ratio=1.50
    dA_pre=.01
    dA_post=.01*ratio 
    trials=10
    
    ### SIMULATION LOOP
    def trial(n): # n is the trial number
        reinit_default_clock()
        clear(True)
    
        eqs_neurons='''
        dv/dt=((ge*(Ee-vr))+Rm*I+(El-v))/taum : volt   
        dge/dt=-ge/taue : 1
        I : amp
        '''
    
        inputs = PoissonGroup(N,rates=lambda t:((.5-.5*cos(2*pi*f*t)))*10*Hz)           
        neurons=NeuronGroup(M,model=eqs_neurons,threshold=vt,reset=vr)
        neurons.I = a*pA
        synapses=Connection(inputs,neurons,'ge',weight=weights)
        neurons.v=vr
    
        S = SpikeMonitor(neurons)
        run(2*second)
        stdp=ExponentialSTDP(synapses,tau_pre,tau_post,dA_pre,-dA_post,wmax=10*weights,interactions='all',update='additive')     
        run(5*second)
        
        phase=zeros((M,200))
        for b in range(0,M):    
            tmp_phase=(S[b]%theta_period)*(360/theta_period)
            phase[b,range(0,len(tmp_phase))] = tmp_phase
            
        return phase
    
    if __name__=='__main__': # This is very important on Windows, otherwise the machine crashes!
        phase = zeros((M,200,trials))
        
        print "This will take approximately 2 minutes."
        pool=multiprocessing.Pool() # uses all available processors b
        results=pool.map(trial,range(trials))
        for i in range(trials):
            phase[:,:,i]=results[i]
        
        ### PLOTTING
        for b in range(0,M):
            m = mean(phase[b,:,:],axis=1)
            st = std(phase[b,:,:],axis=1)/sqrt(trials)
            errorbar(range(0,135), m[range(0,135)], yerr=st[range(0,135)], xerr=None,
                 fmt='-', ecolor=None, elinewidth=None, capsize=3,
                 barsabove=False, lolims=False, uplims=False,
                 xlolims=False, xuplims=False)
        
        title('STDP + Oscillations Simulation')
        xlabel('Spike Number')
        ylabel('Spike Phase (deg)')
        xlim([0, 135])
        ylim([140, 280])
        
        show()