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
|
#!/usr/bin/env python3
"""
Reproduces Figure 2F of
The Use of Hebbian Cell Assemblies for Nonlinear Computation
by Tetzlaff C., Dasgupta S., Kulvicius T. and Wörgötter F.
Sci Rep 5, 12866 (2015).
https://doi.org/10.1038/srep12866
Sebastian Schmitt, 2022
"""
import numpy as np
import matplotlib.pyplot as plt
from brian2 import NeuronGroup, Synapses, StateMonitor, run, defaultclock, ms, second, TimedArray, seed
# random seed that gives curves similar to the ones in the publication
seed(9873487)
# neuron parameters (sigmoidal activation)
beta = 0.03
epsilon = 120
F_max = 100
F_T = 1
tau_u = 1*ms
R = 0.012
# plasticity timescales
tau_ratio = 60
# hebbian
tau_H = 3e4*ms
# synaptic scaling
tau_SS = tau_ratio * tau_H
# synaptic weights
W_max = np.sqrt(tau_ratio*(F_max**2/(F_max - F_T)))
W_ext = W_max
W_input = W_max
W_I = 0.3*W_max
# stimulus
N_units = 100
N_stim_units = 20
stim_A_units_until = N_stim_units
stim_B_units_from = N_units-N_stim_units
# connection probabilities
p_E = 0.1
p_I = 0.2
# paper uses 0.3*ms
DT = 0.5*ms
defaultclock.dt = DT
# duration of a learning trial
lt = 5000*DT
duration = 100*lt
no_input_until = 5*lt
balanced_until = duration/2
# gate balanced presentation of stimulus 1 and 2
balanced = TimedArray([lt_counter*lt < balanced_until for lt_counter in range(int(duration/lt))], dt=lt)
# function used for stimulus (typo in paper, +1 is not part of the argument of sin)
stim_func = TimedArray([100*(np.sin(0.1*(i+1))+1) for i in range(int(duration/DT))], dt=DT)
# gate learning phase of either stimulus 1 or 2
learning_phase = TimedArray([i%10 > 3 for i in range(int(duration/(0.1*lt)))], dt=0.1*lt)
# if not balanced present stimulus A three times more often than stimulus B
stim_A_gate = TimedArray([lt_counter % 2 == 0 if balanced(lt_counter*lt) else lt_counter % 4 in [0,1,2]
for lt_counter in range(int(duration/lt))], dt=lt)
stim_B_gate = TimedArray([lt_counter % 2 == 1 if balanced(lt_counter*lt) else lt_counter % 4 == 3
for lt_counter in range(int(duration/lt))], dt=lt)
# noise is applied also during stimulation
neurons = NeuronGroup(N_units,
"""
F = F_max/(1+exp(beta*(epsilon-u))) : 1
du/dt = (-u + R*(I_E - I_I + W_input*(I_stim_A + I_stim_B)))/tau_u + R*W_ext*20*sqrt((DT/ms)/ms)*xi: 1
I_E : 1
I_I : 1
index : 1 (constant)
stim_units_A = index < stim_A_units_until : boolean
stim_units_B = index >= (stim_B_units_from) : boolean
I_stim_A = learning_phase(t)*int(stim_units_A)*stim_A_gate(t)*stim_func(t) : 1
I_stim_B = learning_phase(t)*int(stim_units_B)*stim_B_gate(t)*stim_func(t) : 1
""",
method = "euler")
neurons.index = range(len(neurons))
# excitatory connections with Hebbian plasticity and synaptic scaling
synapses_E = Synapses(neurons, neurons,
"""
dw/dt = 1/tau_H*F_pre*F_post + 1/tau_SS*(F_T - F_post)*w**2 : 1 (clock-driven)
I_E_post = w*F_pre : 1 (summed)
""",
method="euler"
)
# do not connect between the two populations of stimulated neurons
synapses_E.connect(p=p_E, condition="((j > stim_A_units_until and i >= stim_B_units_from) or (j < stim_B_units_from and i < stim_A_units_until))"
"or ((i > stim_A_units_until and i < stim_B_units_from) and (j > stim_A_units_until and j < stim_B_units_from))")
# fixed weight inhibitory connections
synapses_I = Synapses(neurons, neurons,
"""
w : 1
I_I_post = w*F_pre : 1 (summed)
"""
)
synapses_I.connect(p=p_I)
synapses_I.w = W_I
statemon_neurons = StateMonitor(neurons, ["F", "I_stim_A", "I_stim_B"], record=True, dt=100*defaultclock.dt)
statemon_synapses_E = StateMonitor(synapses_E, "w", record=True, dt=100*defaultclock.dt)
statemon_synapses_for_assembly_analysis = StateMonitor(synapses_E, "w", record=True, dt=lt)
run(duration, report="text")
# threshold saying that synaptic efficacies larger than theta are
# 'strong' and others are 'weak'
theta = 0.5*W_max
in_assembly_A = []
in_assembly_B = []
# traverse through the graph following 'strong' synapses
def go(W, source, units_in_assembly):
units_in_assembly.add(source)
# check all possible targets
for target in range(N_units):
w = W[source][target]
if w > theta:
W[source][target] = 0
go(W, target, units_in_assembly)
# for each learning trial
for ws in statemon_synapses_for_assembly_analysis.w.T:
# construct a full weight matrix
W = np.full((N_units, N_units), np.nan)
W[synapses_E.i[:], synapses_E.j[:]] = ws
for in_assembly, stim_units in zip([in_assembly_A, in_assembly_B],
[range(stim_A_units_until),
range(stim_B_units_from, N_units)]):
units_in_assembly = set()
# start with units that are stimulated
for stim_unit in stim_units:
go(W, stim_unit, units_in_assembly)
in_assembly.append(len(units_in_assembly))
# competitive development of the two competing cell assemblies A and B as a function of the input protocol
fig, ax = plt.subplots()
ax.plot(in_assembly_A, linestyle="None", marker='o', color='orange', label="A")
ax.plot(in_assembly_B, linestyle="None", marker='o', color='olivedrab', label="B")
ax.set_ylim(19, 51)
ax.set_xlim(0, 100)
ax.set_ylabel("Neurons in Cell Assembly [%]")
ax.set_xlabel("Learning Trial")
ax.axvline(balanced_until/lt, linestyle='dashed', color='k')
ax.text(15, 52, " A A", color='orange', fontfamily="monospace", fontsize="xx-large")
ax.text(15, 52, " B B", color='olivedrab', fontfamily="monospace", fontsize="xx-large")
ax.text(65, 52, " 3A 3A", color='orange', fontfamily="monospace", fontsize="xx-large")
ax.text(65, 52, " B B", color='olivedrab', fontfamily="monospace", fontsize="xx-large")
plt.show()
# stimulus, neuronal activity and excitatory weights as function of time
fig, axes = plt.subplots(3, sharex=True)
axes[0].plot(statemon_neurons.I_stim_A[0], label="A", color='orange')
axes[0].plot(statemon_neurons.I_stim_B[-1], label="B", color='olivedrab')
axes[0].legend(loc="upper right")
axes[0].set_title("Stimulus")
axes[1].imshow(statemon_neurons.F, aspect='auto')
axes[1].set_title("Neuron Activity")
axes[1].axhline(stim_A_units_until, linestyle='dashed', color='white')
axes[1].axhline(stim_B_units_from, linestyle='dashed', color='white')
axes[2].imshow(statemon_synapses_E.w, aspect='auto')
axes[2].set_title("Excitatory Weights")
axes[2].set_xticks(range(0, 5000, 250))
axes[2].set_xticklabels(f"{i}" for i in range(0, 100, 5))
axes[2].set_xlabel("Learning Trial")
axes[2].set_xlim(0, 5000)
fig.tight_layout()
plt.show()
|
