Is there a way to define multiple groups of neurons in a loop with all the monitoring functions?

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.