File: examples-frompapers_Diesmann_et_al_1999_longer.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 (243 lines) | stat: -rw-r--r-- 8,817 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
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
.. currentmodule:: brian

.. index::
   pair: example usage; Network
   pair: example usage; Parameters

.. _example-frompapers_Diesmann_et_al_1999_longer:

Example: Diesmann_et_al_1999_longer (frompapers)
================================================

Implementation of synfire chain from Diesmann et al. 1999

Dan Goodman - Dec. 2007

::

    
    #import brian_no_units
    from brian import *
    import time
    
    from brian.library.IF import *
    from brian.library.synapses import *
    
    def minimal_example():
        # Neuron model parameters
        Vr = -70 * mV
        Vt = -55 * mV
        taum = 10 * ms
        taupsp = 0.325 * ms
        weight = 4.86 * mV
        # Neuron model
        equations = Equations('''
            dV/dt = (-(V-Vr)+x)*(1./taum)                            : volt
            dx/dt = (-x+y)*(1./taupsp)                               : volt
            dy/dt = -y*(1./taupsp)+25.27*mV/ms+(39.24*mV/ms**0.5)*xi : volt
            ''')
    
        # Neuron groups
        P = NeuronGroup(N=1000, model=equations,
                      threshold=Vt, reset=Vr, refractory=1 * ms)
    #    P = NeuronGroup(N=1000, model=(dV,dx,dy),init=(0*volt,0*volt,0*volt),
    #                  threshold=Vt,reset=Vr,refractory=1*ms)
    
        Pinput = PulsePacket(t=50 * ms, n=85, sigma=1 * ms)
        # The network structure
        Pgp = [ P.subgroup(100) for i in range(10)]
        C = Connection(P, P, 'y')
        for i in range(9):
            C.connect_full(Pgp[i], Pgp[i + 1], weight)
        Cinput = Connection(Pinput, P, 'y')
        Cinput.connect_full(Pinput, Pgp[0], weight)
        # Record the spikes
        Mgp = [SpikeMonitor(p, record=True) for p in Pgp]
        Minput = SpikeMonitor(Pinput, record=True)
        monitors = [Minput] + Mgp
        # Setup the network, and run it
        P.V = Vr + rand(len(P)) * (Vt - Vr)
        run(100 * ms)
        # Plot result
        raster_plot(showgrouplines=True, *monitors)
        show()
    
    
    # DEFAULT PARAMATERS FOR SYNFIRE CHAIN
    # Approximates those in Diesman et al. 1999
    model_params = Parameters(
        # Simulation parameters
        dt=0.1 * ms,
        duration=100 * ms,
        # Neuron model parameters
        taum=10 * ms,
        taupsp=0.325 * ms,
        Vt= -55 * mV,
        Vr= -70 * mV,
        abs_refrac=1 * ms,
        we=34.7143,
        wi= -34.7143,
        psp_peak=0.14 * mV,
        # Noise parameters
        noise_neurons=20000,
        noise_exc=0.88,
        noise_inh=0.12,
        noise_exc_rate=2 * Hz,
        noise_inh_rate=12.5 * Hz,
        computed_model_parameters="""
        noise_mu = noise_neurons * (noise_exc * noise_exc_rate - noise_inh * noise_inh_rate ) * psp_peak * we
        noise_sigma = (noise_neurons * (noise_exc * noise_exc_rate + noise_inh * noise_inh_rate ))**.5 * psp_peak * we
        """
        )
    
    # MODEL FOR SYNFIRE CHAIN
    # Excitatory PSPs only
    def Model(p):
        equations = Equations('''
            dV/dt = (-(V-p.Vr)+x)*(1./p.taum)                          : volt
            dx/dt = (-x+y)*(1./p.taupsp)                               : volt
            dy/dt = -y*(1./p.taupsp)+25.27*mV/ms+(39.24*mV/ms**0.5)*xi : volt
            ''')
        return Parameters(model=equations, threshold=p.Vt, reset=p.Vr, refractory=p.abs_refrac)
    
    default_params = Parameters(
        # Network parameters
        num_layers=10,
        neurons_per_layer=100,
        neurons_in_input_layer=100,
        # Initiating burst parameters
        initial_burst_t=50 * ms,
        initial_burst_a=85,
        initial_burst_sigma=1 * ms,
        # these values are recomputed whenever another value changes
        computed_network_parameters="""
        total_neurons = neurons_per_layer * num_layers
        """,
        # plus we also use the default model parameters
        ** model_params
        )
    
    # DEFAULT NETWORK STRUCTURE
    # Single input layer, multiple chained layers
    class DefaultNetwork(Network):
        def __init__(self, p):
            # define groups
            chaingroup = NeuronGroup(N=p.total_neurons, **Model(p))
            inputgroup = PulsePacket(p.initial_burst_t, p.neurons_in_input_layer, p.initial_burst_sigma)
            layer = [ chaingroup.subgroup(p.neurons_per_layer) for i in range(p.num_layers) ]
            # connections
            chainconnect = Connection(chaingroup, chaingroup, 2)
            for i in range(p.num_layers - 1):
                chainconnect.connect_full(layer[i], layer[i + 1], p.psp_peak * p.we)
            inputconnect = Connection(inputgroup, chaingroup, 2)
            inputconnect.connect_full(inputgroup, layer[0], p.psp_peak * p.we)
            # monitors
            chainmon = [SpikeMonitor(g, True) for g in layer]
            inputmon = SpikeMonitor(inputgroup, True)
            mon = [inputmon] + chainmon
            # network
            Network.__init__(self, chaingroup, inputgroup, chainconnect, inputconnect, mon)
            # add additional attributes to self
            self.mon = mon
            self.inputgroup = inputgroup
            self.chaingroup = chaingroup
            self.layer = layer
            self.params = p
    
        def prepare(self):
            Network.prepare(self)
            self.reinit()
    
        def reinit(self, p=None):
            Network.reinit(self)
            q = self.params
            if p is None: p = q
            self.inputgroup.generate(p.initial_burst_t, p.initial_burst_a, p.initial_burst_sigma)
            self.chaingroup.V = q.Vr + rand(len(self.chaingroup)) * (q.Vt - q.Vr)
    
        def run(self):
            Network.run(self, self.params.duration)
    
        def plot(self):
            raster_plot(ylabel="Layer", title="Synfire chain raster plot",
                       color=(1, 0, 0), markersize=3,
                       showgrouplines=True, spacebetweengroups=0.2, grouplinecol=(0.5, 0.5, 0.5),
                       *self.mon)
    
    def estimate_params(mon, time_est):
        # Quick and dirty algorithm for the moment, for a more decent algorithm
        # use leastsq algorithm from scipy.optimize.minpack to fit const+Gaussian
        # http://www.scipy.org/doc/api_docs/SciPy.optimize.minpack.html#leastsq
        i, times = zip(*mon.spikes)
        times = array(times)
        times = times[abs(times - time_est) < 15 * ms]
        if len(times) == 0:
            return (0, 0 * ms)
        better_time_est = times.mean()
        times = times[abs(times - time_est) < 5 * ms]
        if len(times) == 0:
            return (0, 0 * ms)
        return (len(times), times.std()*second)
    
    def single_sfc():
        net = DefaultNetwork(default_params)
        net.run()
        net.plot()
    
    def state_space(grid, neuron_multiply, verbose=True):
        amin = 0
        amax = 100
        sigmamin = 0. * ms
        sigmamax = 3. * ms
    
        params = default_params()
        params.num_layers = 1
        params.neurons_per_layer = params.neurons_per_layer * neuron_multiply
    
        net = DefaultNetwork(params)
    
        i = 0
        # uncomment these 2 lines for TeX labels
        #import pylab
        #pylab.rc_params.update({'text.usetex': True}) 
        if verbose:
            print "Completed:"
        start_time = time.time()
        figure()
        for ai in range(grid + 1):
            for sigmai in range(grid + 1):
                a = int(amin + (ai * (amax - amin)) / grid)
                if a > amax: a = amax
                sigma = sigmamin + sigmai * (sigmamax - sigmamin) / grid
                params.initial_burst_a, params.initial_burst_sigma = a, sigma
                net.reinit(params)
                net.run()
                (newa, newsigma) = estimate_params(net.mon[-1], params.initial_burst_t)
                newa = float(newa) / float(neuron_multiply)
                col = (float(ai) / float(grid), float(sigmai) / float(grid), 0.5)
                plot([sigma / ms, newsigma / ms], [a, newa], color=col)
                plot([sigma / ms], [a], marker='.', color=col, markersize=15)
                i += 1
                if verbose:
                    print str(int(100. * float(i) / float((grid + 1) ** 2))) + "%",
            if verbose:
                print
        if verbose:
            print "Evaluation time:", time.time() - start_time, "seconds"
        xlabel(r'$\sigma$ (ms)')
        ylabel('a')
        title('Synfire chain state space')
        axis([sigmamin / ms, sigmamax / ms, amin, amax])
    
    
    minimal_example()
    #print 'Computing SFC with multiple layers'
    #single_sfc()
    #print 'Plotting SFC state space'
    #state_space(3,1)
    #state_space(8,10)
    #state_space(10,50)
    #state_space(10,150)
    #show()