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
|
.. currentmodule:: brian
.. index::
pair: example usage; load
pair: example usage; plot
pair: example usage; run
pair: example usage; log
pair: example usage; show
pair: example usage; network_operation
pair: example usage; sum
pair: example usage; semilogy
pair: example usage; argsort
pair: example usage; Connection
pair: example usage; subplot
pair: example usage; SpikeMonitor
pair: example usage; linspace
pair: example usage; exp
pair: example usage; NeuronGroup
pair: example usage; nonzero
.. _example-frompapers-computing with neural synchrony-olfaction_Fig11B_olfaction_stdp_testing:
Example: Fig11B_olfaction_stdp_testing (frompapers/computing with neural synchrony/olfaction)
=============================================================================================
Brette R (2012). Computing with neural synchrony. PLoS Comp Biol. 8(6): e1002561. doi:10.1371/journal.pcbi.1002561
------------------------------------------------------------------------------------------------------------------
Figure 11B. Learning to recognize odors.
Caption (Fig. 11B). After learning, responses
of postsynaptic neurons, ordered by tuning ratio, to odor A (blue) and odor B (red),
with an increasing concentration (0.1 to 10, where 1 is odor
concentration in the learning phase).
Run the other file first: Fig11B_olfaction_stdp_learning.py
::
from brian import *
from params import *
import numpy
from scipy.sparse import lil_matrix
bmin,bmax=-7,-1
# Loads information from the STDP simulation
t,odor=numpy.load("stimuli.npy").T
W=numpy.load("weights.npy")
spikes_out=numpy.load("spikesout.npy")
weights=W[-1,:,:] # final weights
# Analyze selectivity at the end of the STDP simulation
ispikes=spikes_out[:,0] # indexes of neurons that spiked
tspikes=spikes_out[:,1] # spike timings
# Select only the end of the STDP simulation
end=tspikes>.8*max(tspikes)
ispikes=ispikes[end]
tspikes=tspikes[end]
odors=odor[digitize(tspikes,t)-1] # odor (0/1) presented at the time of spikes
tuning=zeros(30) # Tuning ratio of the postsynaptic neurons
n0,n1=zeros(30),zeros(30) # number of spikes for odor 0 and for odor 1
for k in range(len(tuning)):
o=odors[ispikes==k]
n0[k]=sum(o==0)
n1[k]=sum(o==1)
tuning[k]=n0[k]*1./(n0[k]+n1[k])
# Sort the postsynaptic neurons by odor tuning
weights=weights[:,argsort(tuning)]
'''
Run the simulation
'''
def odor(N):
# Returns a random vector of binding constants
return 10**(rand(N)*(bmax-bmin)+bmin)
def hill_function(c,K=1.,n=3.):
'''
Hill function:
* c = concentration
* K = half activation constant (choose K=1 for relative concentrations)
* n = Hill coefficient
'''
return (c**n)/(c**n+K**n)
N=5000 # number of receptors
seed(31415) # Get the same neurons every time
intensity=3000.
# Odor plumes
tau_plume=75*ms
eq_plumes='''
dx/dt=-x/tau_plume+(2./tau_plume)**.5*xi : 1
y=clip(x,0,inf) : 1
'''
plume=NeuronGroup(2,model=eq_plumes) # 1 odor
# Receptor neurons
Fmax=40*Hz # maximum firing rate
tau=20*ms
Imax=1/(1-exp(-1/(Fmax*tau))) # maximum input current
eq_receptors='''
dv/dt=(Imax*hill_function(c)-v)/tau : 1
c : 1 # concentrations (relative to activation constant)
'''
receptors=NeuronGroup(N,model=eq_receptors,threshold=1,reset=0)
@network_operation
def odor_to_nose():
# Send odor plume to the receptors
receptors.c=I1*c1*clip(plume.x[0],0,Inf)+I2*c2*clip(plume.x[0],0,Inf)
odors=[odor(N),odor(N)]
c1,c2=odors
# Decoder neurons
M=len(tuning)
eq_decoders='''
dv/dt=-v/taud + sigma*(2/taud)**.5*xi : 1
'''
decoders=NeuronGroup(M,model=eq_decoders,threshold=1,reset=0)
S2=SpikeMonitor(decoders)
# Synapses
syn=Connection(receptors,decoders,'v')
for i in range(len(decoders)):
for j in weights[:,i].nonzero()[0]:
syn[j,i]=weights[j,i]
# Run
I1,I2=intensity,0
print "Started"
# Odor A, increasing concentration
for I1 in intensity*exp(linspace(log(.1),log(10),20)):
run(1*second,report="text")
I1=0
# Odor B, increasing concentration
for I2 in intensity*exp(linspace(log(.1),log(10),20)):
run(1*second,report="text")
# Figure (11B)
spikes=array(S2.spikes) # i,t
n,t=spikes[:,0],spikes[:,1]
subplot(211) # Raster plot
plot(t,n,'k.')
subplot(212) # Odor concentrations
semilogy(linspace(0,20,20),exp(linspace(log(.1),log(10),20)),'b')
semilogy(linspace(20,40,20),exp(linspace(log(.1),log(10),20)),'r')
plot([0,40],[1,1],'k--')
show()
|
