Hi. I’m trying to define many groups of neurons and need to define them in a loop. My problem is that when I append them in a list, it seems that I cannot monitor the variables that I want. I would really appreciate it if someone could help me in this regard. Thank you.
Hi LBENE,
i often use python’s list comprehension to solve this:
voltage_monitors = [StateMonitor(pop,'v',record=True) for pop in pop_list]
spike_monitors = [SpikeMonitor(pop) for pop in pop_list]
rate_monitors = [PopulationRateMonitor(pop) for pop in pop_list]
where pop_list is just a list of neuron groups
1 Like
Thank you very much adam. Are the synapses between each two neuron groups defined the same way?
I haven’t tried monitoring synapses yet, but yeah I think if you collected synapses like:
a_syn = Synapses(...)
b_syn = Synapses(...)
c_syn = Synapses(...)
synapse_list = [a_syn, b_syn, c_syn]
or
synapse_list = []
for i in range(3):
synapse_list.append( Synapses(...) )
then you could do
synpase_w_monitors = [StateMonitor(S, 'w', record=True) for S in synapse_list]
1 Like
I tried expanding the code to a larger network and the same thing happened again. StateMonitor records nothing. This is the related part of the code (e: excitatory, ffi: feedforward inhibitory, fbi: feedfback inhibitory):
...
G_e = []
G_ffi = []
G_fbi = []
S_e_e = []
S_e_ffi = []
S_e_fbi = []
S_ffi_e = []
S_ffi_ffi = []
S_fbi_e = []
S_fbi_fbi = []
for rc1 in range(number_of_regions):
...
G_e.append(NeuronGroup(net1_attr[rc1].N_e, eqs_e, threshold="v > 0*mV", refractory="v > 0*mV", method="exponential_euler"))
...
G_ffi.append(NeuronGroup(net1_attr[rc1].N_ffi, eqs_ffi, threshold="v > 0*mV", refractory="v > 0*mV", method="exponential_euler"))
G_fbi.append(NeuronGroup(net1_attr[rc1].N_fbi, eqs_fbi, threshold="v > 0*mV", refractory="v > 0*mV", method="exponential_euler"))
S_e_e.append(Synapses(G_e[rc1], G_e[rc1], "w : volt", on_pre="v_post += w"))
S_e_ffi.append(Synapses(G_e[rc1], G_ffi[rc1], "w : volt", on_pre="v_post += w"))
S_e_fbi.append(Synapses(G_e[rc1], G_fbi[rc1], "w : volt", on_pre="v_post += w"))
S_ffi_e.append(Synapses(G_ffi[rc1], G_e[rc1], "w : volt", on_pre="v_post -= w"))
S_ffi_ffi.append(Synapses(G_ffi[rc1], G_ffi[rc1], "w : volt", on_pre="v_post -= w"))
S_fbi_e.append(Synapses(G_fbi[rc1], G_e[rc1], "w : volt", on_pre="v_post -= w"))
S_fbi_fbi.append(Synapses(G_fbi[rc1], G_fbi[rc1], "w : volt", on_pre="v_post -= w"))
...
S_e_e[rc1].connect(i=i1, j=post1)
...
S_e_ffi[rc1].connect(i=i1, j=post1)
...
S_e_fbi[rc1].connect(i=i1, j=post1)
...
S_e_e[rc1].w = scale1*np.random.normal(PP_m, STD1)*mV
S_e_ffi[rc1].w = scale1*np.random.normal(PPV_m, STD1)*mV
S_e_fbi[rc1].w = scale1*np.random.normal(PS_m, STD1)*mV
S_e_e[rc1].delay = delay_e_e*ms
S_e_ffi[rc1].delay = delay_e_ffi*ms
S_e_fbi[rc1].delay = delay_e_fbi*ms
>>>same approach for ffi and fbi<<<
...
S_e_e_o = []
S_e_ffi_o = []
S_ffi_e_o = []
S_ffi_ffi_o = []
for rc1 in range(number_of_regions):
S_e_e_o.append([])
S_e_ffi_o.append([])
S_ffi_e_o.append([])
S_ffi_ffi_o.append([])
for rc2 in range(number_of_regions):
if rc1 != rc2:
S_e_e_o[rc1].append(Synapses(G_e[rc1], G_e[rc2], "w : volt", on_pre="v_post += w"))
S_e_ffi_o[rc1].append(Synapses(G_e[rc1], G_ffi[rc2], "w : volt", on_pre="v_post += w"))
S_ffi_e_o[rc1].append(Synapses(G_ffi[rc1], G_e[rc2], "w : volt", on_pre="v_post -= w"))
S_ffi_ffi_o[rc1].append(Synapses(G_ffi[rc1], G_ffi[rc2], "w : volt", on_pre="v_post -= w"))
S_e_e_o[rc1][rc2].connect(p=prob_t[0])
S_e_ffi_o[rc1][rc2].connect(p=prob_t[1])
S_ffi_e_o[rc1][rc2].connect(p=prob_t[2])
S_ffi_ffi_o[rc1][rc2].connect(p=prob_t[3])
S_e_e_o[rc1][rc2].w = scale1*np.random.normal(PP_m, STD1)*mV
S_e_ffi_o[rc1][rc2].w = scale1*np.random.normal(PPV_m, STD1)*mV
S_ffi_e_o[rc1][rc2].w = scale1*np.random.normal(PVP_m, STD1)*mV
S_ffi_ffi_o[rc1][rc2].w = scale1*np.random.normal(PVPV_m, STD1)*mV
S_e_e_o[rc1][rc2].delay = 20.*ms
S_e_ffi_o[rc1][rc2].delay = 13.*ms
S_ffi_e_o[rc1][rc2].delay = 17.*ms
S_ffi_ffi_o[rc1][rc2].delay = 17.*ms
else:
S_e_e_o[rc1].append([])
S_e_ffi_o[rc1].append([])
S_ffi_e_o[rc1].append([])
S_ffi_ffi_o[rc1].append([])
M_v_e = [StateMonitor(g_e,'v',record=True) for g_e in G_e]
M_v_ffi = [StateMonitor(g_ffi,'v',record=True) for g_ffi in G_ffi]
M_v_fbi = [StateMonitor(g_fbi,'v',record=True) for g_fbi in G_fbi]
run(60*second)
Hi @LBENE , please see the remarks in the documentation here: Running a simulation — Brian 2 2.5.0.1 documentation If your objects (e.g. monitors) are in a list, you have to manually add them to a Network, a simple run will not include them. See also my reply here: Unable to create proper synapses connections using a loop - #2 by mstimberg
2 Likes
Oh, and one more thing: it is much more efficient to create a minimal number of NeuronGroup and Synapses objects. In particular, this will drastically reduce the preparation/compilation time. From your code example, it seems that all the neurons stored in G_e share the same equations, they could therefore all go into a single NeuronGroup (same for G_ffi and G_fbi). You can always create subgroups to divide them up into smaller units. Similarly, you could then have considerable fewer Synapses objects, since you only need one per type of connection.
A technique that is not widely used but can lead to very readable code is to make use of “labels” in the neuron description, e.g. add region : integer (constant) to the equations. Then, assuming you have a single NeuronGroup G_e storing all layers of excitatory neurons, you’ll set them with something like (inspired from the code you posted):
region_sizes = [net1_attr[rc1].N_e for rc1 in range(number_regions)]
indices = np.cumsum(np.hstack([0, region_sizes]))
for r in range(number_of_regions):
G_e.region[indices[r]:indices[r+1]] = r
You can then refer to the region of cells when you create synapses. For example, you could create connections between the excitatory neurons in each region with a single connect command:
S_e_e = Synapses(G_e, G_e, ...)
S_e_e.connect("region_pre == region_post", p=0.1)
The only downside of this approach is that it is not the most efficient way to connect the synapses: it will check for all possible pairs of neurons whether they are in the same region. If you have a very big network, this will be quite slow.
1 Like
Wow! Thank you very much for all the insights. The network is not that big right now: five regions with 250 neurons each. The number of neurons in each region will be increased later on to something like 1000.
Also, I previously overlooked the part about subgroups. After this, I wanted to define like 75 groups for each region representing 25 cortical columns. This subgroup really avoids going on that path.
1 Like
Hi again. I tried breaking the groups in smaller pieces and it apparently got very slow. I tried printing some variables and it seems for some reason synapses are defined very slow. This is the related part of the code:
for rc1 in range(number_of_regions):
...
for i1 in range(20):
G1 = G_py[rc1][14*i1:14*(i1+1)]
for i2 in range(20):
G2 = G_py[rc1][14*i2:14*(i2+1)]
syn_temp = Synapses(G1, G2, "w : volt", on_pre="v_post += w")
if i1 == i2:
syn_temp.connect(p=.65)
syn_temp.w = scale1*np.random.normal(PYPY_m, STD1)*mV
syn_temp.delay = synaptic_time_constant_py_py*ms
else:
syn_temp.connect(p=.1)
syn_temp.w = scale1*np.random.normal(PYPY_m, STD1)*mV
syn_temp.delay = (neurons_distance(i1, i2) / conduction_velocity_py_py + synaptic_time_constant_py_py)*ms
S_all_intra.append(syn_temp)
if i2%2 == 0:
G3 = G_pv[rc1][5*int(np.floor(i2/2)):5*int(np.floor(i2/2))+3]
else:
G3 = G_pv[rc1][5*int(np.floor(i2/2))+3:5*int(np.floor(i2/2))+5]
syn_temp = Synapses(G1, G3, "w : volt", on_pre="v_post += w")
if i1 == i2:
syn_temp.connect(p=.65)
syn_temp.w = scale1*np.random.normal(PYPV_m, STD1)*mV
syn_temp.delay = synaptic_time_constant_py_pv*ms
else:
syn_temp.connect(p=.1)
syn_temp.w = scale1*np.random.normal(PYPV_m, STD1)*mV
syn_temp.delay = (neurons_distance(i1, i2) / conduction_velocity_py_pv + synaptic_time_constant_py_pv)*ms
S_all_intra.append(syn_temp)
G4 = G_sst[rc1][i2]
syn_temp = Synapses(G1, G4, "w : volt", on_pre="v_post += w")
if i1 == i2:
syn_temp.connect(p=.65)
syn_temp.w = scale1*np.random.normal(PYSST_m, STD1)*mV
syn_temp.delay = synaptic_time_constant_py_sst*ms
else:
syn_temp.connect(p=.1)
syn_temp.w = scale1*np.random.normal(PYSST_m, STD1)*mV
syn_temp.delay = (neurons_distance(i1, i2) / conduction_velocity_py_sst + synaptic_time_constant_py_sst)*ms
S_all_intra.append(syn_temp)
if i1%2 == 0:
G1 = G_pv[rc1][5*int(np.floor(i1/2)):5*int(np.floor(i1/2))+3]
else:
G1 = G_pv[rc1][5*int(np.floor(i1/2))+3:5*int(np.floor(i1/2))+5]
for i2 in range(20):
G2 = G_py[rc1][14*i2:14*(i2+1)]
syn_temp = Synapses(G1, G2, "w : volt", on_pre="v_post -= w")
if i1 == i2:
syn_temp.connect(p=.62)
syn_temp.w = scale1*np.random.normal(PVPY_m, STD1)*mV
syn_temp.delay = synaptic_time_constant_pv_py*ms
else:
syn_temp.connect(p=.07)
syn_temp.w = scale1*np.random.normal(PVPY_m, STD1)*mV
syn_temp.delay = (neurons_distance(i1, i2) / conduction_velocity_pv_py + synaptic_time_constant_pv_py)*ms
S_all_intra.append(syn_temp)
if i2%2 == 0:
G3 = G_pv[rc1][5*int(np.floor(i2/2)):5*int(np.floor(i2/2))+3]
else:
G3 = G_pv[rc1][5*int(np.floor(i2/2))+3:5*int(np.floor(i2/2))+5]
syn_temp = Synapses(G1, G3, "w : volt", on_pre="v_post -= w")
if i1 == i2:
syn_temp.connect(p=.8)
syn_temp.w = scale1*np.random.normal(PVPV_m, STD1)*mV
syn_temp.delay = synaptic_time_constant_pv_pv*ms
else:
syn_temp.connect(p=.07)
syn_temp.w = scale1*np.random.normal(PVPV_m, STD1)*mV
syn_temp.delay = (neurons_distance(i1, i2) / conduction_velocity_pv_pv + synaptic_time_constant_pv_pv)*ms
S_all_intra.append(syn_temp)
G4 = G_sst[rc1][i2]
syn_temp = Synapses(G1, G4, "w : volt", on_pre="v_post -= w")
if i1 == i2:
syn_temp.connect(p=.8)
syn_temp.w = scale1*np.random.normal(PVSST_m, STD1)*mV
syn_temp.delay = synaptic_time_constant_pv_sst*ms
else:
syn_temp.connect(p=.07)
syn_temp.w = scale1*np.random.normal(PVSST_m, STD1)*mV
syn_temp.delay = (neurons_distance(i1, i2) / conduction_velocity_pv_sst + synaptic_time_constant_pv_sst)*ms
S_all_intra.append(syn_temp)
G1 = G_sst[rc1][i1]
for i2 in range(20):
G2 = G_py[rc1][14*i2:14*(i2+1)]
syn_temp = Synapses(G1, G2, "w : volt", on_pre="v_post -= w")
if i1 == i2:
syn_temp.connect(p=.71)
syn_temp.w = scale1*np.random.normal(SSTPY_m, STD1)*mV
syn_temp.delay = synaptic_time_constant_sst_py*ms
else:
syn_temp.connect(p=0.)
syn_temp.w = scale1*np.random.normal(SSTPY_m, STD1)*mV
syn_temp.delay = (neurons_distance(i1, i2) / conduction_velocity_sst_py + synaptic_time_constant_sst_py)*ms
S_all_intra.append(syn_temp)
if i2%2 == 0:
G3 = G_pv[rc1][5*int(np.floor(i2/2)):5*int(np.floor(i2/2))+3]
else:
G3 = G_pv[rc1][5*int(np.floor(i2/2))+3:5*int(np.floor(i2/2))+5]
syn_temp = Synapses(G1, G3, "w : volt", on_pre="v_post -= w")
if i1 == i2:
syn_temp.connect(p=.8)
syn_temp.w = scale1*np.random.normal(SSTPV_m, STD1)*mV
syn_temp.delay = synaptic_time_constant_sst_pv*ms
else:
syn_temp.connect(p=0.)
syn_temp.w = scale1*np.random.normal(SSTPV_m, STD1)*mV
syn_temp.delay = (neurons_distance(i1, i2) / conduction_velocity_sst_pv + synaptic_time_constant_sst_pv)*ms
S_all_intra.append(syn_temp)
I would really appreciate it if you could guide me in this regard.
I rewrote the code like:
for rc1 in range(number_of_regions):
S_py_py.append(Synapses(G_py[rc1], G_py[rc1], "w : volt", on_pre="v_post += w"))
...
for i1 in range(20):
g1_0 = range(14*i1,14*(i1+1))
for i2 in range(20):
g2_0 = range(14*i2,14*(i2+1))
G1_0 = np.repeat(g1_0, np.size(g2_0))
G2_0 = np.tile(g2_0, np.size(g1_0))
if i1 == i2:
S_py_py[rc1].connect(i=G1_0,j=G2_0,p=.65)
S_py_py[rc1].w[G1_0,G2_0] = scale1*np.random.normal(PYPY_m, STD1)*mV
S_py_py[rc1].delay[G1_0,G2_0] = synaptic_time_constant_py_py*ms
else:
S_py_py[rc1].connect(i=G1_0,j=G2_0,p=.1)
S_py_py[rc1].w[G1_0,G2_0] = scale1*np.random.normal(PYPY_m, STD1)*mV
S_py_py[rc1].delay[G1_0,G2_0] = (neurons_distance(i1, i2) / conduction_velocity_py_py + synaptic_time_constant_py_py)*ms
And it seems fast now. My only problem is that I don’t know if the way I defined connections, weights, and delays works fine. In other words, can this setup be used for having various weights and delays inside a synaptic object?
Hi. Thanks for the example code, that is really helpful. I’ll have a closer look, but I think whether this works correctly depends on whether the S_py_py[rc1].w = ... and ...delay = ... lines are setting single values, or an array of values. In other words should all connections within a column, or all connections between two specific columns, share the same random weight? This seems to be the case in your current code (assuming that PYPY_m and STD1 are scalar values which means that np.random.normal returns a single random value).
1 Like