-
Notifications
You must be signed in to change notification settings - Fork 1
/
GutigSompolinsky2009.py
153 lines (134 loc) · 5.2 KB
/
GutigSompolinsky2009.py
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
from brian import *
class TimeWarpModel(object):
'''
A simple neuron model for testing the "time-warp invariance" with
conductance based or current based synapses. The neuron receives balanced
excitatory and inhibitory input from a random spike train. The same spike
train can be fed into the model with different time warps.
'''
def __init__(self, conductance_based=True):
'''
Create a new model with conductance based or current based synapses
'''
# Model parameters
E_e = 5
E_i = -1
E_L = 0
g_L = 1/(100.0*msecond)
tau_syn = 1*ms
N_ex = 250
N_inh = 250
self.N = N_ex + N_inh
# Equations
if conductance_based:
eqs = '''
dV/dt = -(V - E_L) * g_L - I_syn : 1
I_syn = I_ge + I_gi : second**-1
I_ge = (V - E_e) * g_e : second**-1
I_gi = (V - E_i) * g_i : second**-1
dg_e/dt = -g_e/tau_syn : second**-1
dg_i/dt = -g_i/tau_syn : second**-1
'''
else:
eqs = '''
dV/dt = -(V - E_L) * g_L - I_syn : 1
I_syn = -5 * g_e + g_i : second**-1
dg_e/dt = -g_e/tau_syn : second**-1
dg_i/dt = -g_i/tau_syn : second**-1
'''
# for simpler voltage traces: no spiking
neuron = NeuronGroup(1, model=eqs, threshold=None)
# every input neuron fires once in a random interval
self.unwarped_spiketimes = [(i, t * 250 * ms) for i, t in
zip(range(0, self.N), rand(self.N))]
# final spiketimes will be set in the run function
self.input = SpikeGeneratorGroup(self.N, [])
e_input = self.input.subgroup(N_ex)
i_input = self.input.subgroup(N_inh)
e_conn = Connection(e_input, neuron, 'g_e',
weight=6 / (N_ex * tau_syn))
i_conn = Connection(i_input, neuron, 'g_i',
weight=5 * 6 / (N_ex * tau_syn))
# record membrane potential
self.monitor = StateMonitor(neuron, varname='V', record=True)
# putting everything together
self.net = Network(neuron, self.input, e_conn, i_conn, self.monitor)
def run(self, beta=1.0):
'''
Run the network with the original spike train warped by a certain factor
beta. Beta > 1 corresponds to an extended and beta < 1 to a shrinked
input spike train.
'''
self.net.reinit()
#warp spike train in time
self.input.spiketimes = [(i, beta*t)
for i, t in self.unwarped_spiketimes]
self.net.run(beta * 250*ms)
#Return the voltage trace
return (self.monitor.times, self.monitor[0])
if __name__ == '__main__':
cond_model = TimeWarpModel(True)
curr_model = TimeWarpModel(False)
N = cond_model.N
# #########################################################################
# Reproduce Fig. 2 from Gutig and Sompolinsky (2009)
# #########################################################################
beta = 2.0
times1, v1 = cond_model.run(beta=1.0)
times2, v2 = cond_model.run(beta=beta)
maxtime = 250 * beta
subplot(4, 1, 1)
(neurons, times) = zip(*cond_model.unwarped_spiketimes)
plot(array(times) / ms, neurons, 'g.')
axis([0, maxtime, 0, N])
xticks([])
yticks([])
title('Time-warp-invariant voltage traces (conductance-based)')
subplot(4, 1, 2)
plot(times1 / ms, v1, 'g')
axis([0, maxtime, -1.5, 1.5])
xticks([])
yticks([])
subplot(4, 1, 3)
plot(array(times) * beta / ms, neurons, 'b.')
axis([0, maxtime, 0, N])
xticks([])
yticks([1, 500])
subplot(4, 1, 4)
plot(times2 / ms, v2, 'b')
plot(times1 / ms * beta, v1, 'g')
axis([0, maxtime, -1.5, 1.5])
xlabel('Time (ms)')
xticks([0, 250, 500])
yticks([-1, 1])
show()
# #########################################################################
# Reproduce Fig. 3(C) from Gutig and Sompolinsky (2009), but for random
# spike trains and not in a speech recognition task
# #########################################################################
betas = arange(0.2, 3.1, 0.1)
#betas = array([1.0, 2.0])
cond_results = []
curr_results = []
for beta in betas:
print 'Testing warp factor %.1f' % beta
cond_results.append(cond_model.run(beta))
curr_results.append(curr_model.run(beta))
figure()
colors = mpl.cm.gist_earth((betas - betas[0]) / (betas[-1] - betas[0]))
lookup = dict(zip(betas, colors))
for beta, cond_result, curr_result in zip(betas, cond_results,
curr_results):
times_cond, v_cond = cond_result
times_curr, v_curr = curr_result
subplot(1,2,1)
plot(times_cond / ms / beta, v_cond, color=lookup[beta])
axis([0, 250, -1.5, 1.5])
subplot(1,2,2)
plot(times_curr / ms / beta, v_curr, color=lookup[beta])
axis([0, 250, -1.5, 1.5])
subplot(1,2,1)
title('conductance based')
subplot(1,2,2)
title('current based')
show()