External Voltage/Some Neurons not spiking

Support
1

Description of problem

Hope you are doing well. Im trying to give an external voltage to a random subset of neurons, and trying to see how the signal propagates. I have given the external voltage to a random 10% subset, and seems like that subset spikes as expected. But I have connected that Inhibitory group to another excitatory group, but the excitatory group isn’t spiking and nor is the inhibitory group crossing the threshold voltage.
Minimal code to reproduce problem


net = Network()

num_excitatory = 103
num_inhibitory = 103
num_dopamine = 103
pic_dur = 4*ms
_start = 0
_stop = 1_000 * ms
_step = 0.2 * ms
times = np.linspace(_start, _stop, int((_stop - _start) / _step) + 1)
values = 100 * volt * np.ones(len(times))
V_ext = TimedArray(values, dt=_step)
eqs = '''dv/dt = (ge * (Ee - v) + El - v+(int(receives_input)*V_ext(t))) / taum : volt (unless refractory)
           dge/dt = -ge / taue : 1
           receives_input: boolean (constant)'''

eqs_S_E = '''mode: 1
                        w:1
                         dc/dt = -c / tauc : 1 (clock-driven)
                         dd/dt = -d / taud : 1 (clock-driven)
                         ds/dt = mode * c * d / taus : 1 (clock-driven)
                         dApre/dt = -Apre / taupre : 1 (event-driven)
                         dApost/dt = -Apost / taupost : 1 (event-driven)'''

eqs_pre = '''ge += (s+w)
                          Apre += dApre
                          c = clip(c + mode * Apost, -gmax, gmax)
                          s = clip(s + (1-mode) * Apost, -gmax, gmax)
                          '''

eqs_post = '''Apost += dApost_E
    c = clip(c + mode * Apre, -Egmax, Egmax)
    s = clip(s + (1-mode) * Apre, 0, Egmax)
    '''
eqs_S_I = '''mode: 1
    w:1
    dc/dt = -c / tauc : 1 (clock-driven)
    dd/dt = -d / taud : 1 (clock-driven)
    ds/dt = mode * c * d / taus : 1 (clock-driven)
    dApre/dt = - (Apre + dApre_I/15) / taupre : 1 (event-driven)
    dApost/dt = - (Apost + dApre_I/15) / taupre : 1 (event-driven)'''

eqs_pre_I = '''
    ge -= (s+w)
    Apre += dApre_I
    c = clip(c + mode * Apost, -Igmax, Igmax)
    s = clip(s + (1-mode) * Apost, 0, Igmax)
    '''

eqs_post_I = '''Apost += dApre_I
    c = clip(c + mode * Apre, -Igmax, Igmax)
    s = clip(s + (1-mode) * Apre, 0, Igmax)
    '''

Excitatory  =  NeuronGroup(num_excitatory, eqs, threshold='v>vt', reset='v = vr',  method='euler',refractory= 2 * ms)
Inhibitory = NeuronGroup(num_inhibitory, eqs, threshold='v>vt', reset='v = vr',  method='euler',refractory= 2 * ms)
dopamine = NeuronGroup(num_dopamine, '''dv/dt =  - v / (100 * ms) : volt (unless refractory)''',threshold='v>100 * mV', refractory= 2 * ms,reset='v=0 * mV', method='linear')
Excitatory.v = vr
Inhibitory.v = vr
selected_neurons = np.random.choice(num_neurons, int(0.1*num_neurons), replace=False)
for neuron_index in selected_neurons:
    print(neuron_index)
    Inhibitory[neuron_index].receives_input = True
dopamine.v = vr
net.add(Excitatory)
net.add(Inhibitory)
net.add(dopamine)

#synapse = Synapses(spike_input,Excitatory ,model='''s: volt''',on_pre='v += s')
#synapse.connect(i=[0, 1], j=[0, 1])
#synapse.s = 100. * mV
#net.add(synapse)
C_E_E = Synapses(Excitatory,Excitatory , model=eqs_S_E, on_pre=eqs_pre, on_post=eqs_post, method='euler')
C_E_I = Synapses(Excitatory ,Inhibitory ,model=eqs_S_E, on_pre=eqs_pre, on_post=eqs_post, method='euler')
C_I_I = Synapses(Inhibitory, Inhibitory, model=eqs_S_I, on_pre=eqs_pre_I, on_post=eqs_post_I, method='euler')
C_I_E = Synapses(Inhibitory, Excitatory, model=eqs_S_I, on_pre=eqs_pre_I, on_post=eqs_post_I, method='euler')
C_D_I = Synapses(dopamine, Inhibitory,model='''w:1''',on_pre='''d_post += epsilon_dopa * (dopamine_times(t) * 2 - 1)''',method='exact')
C_D_E = Synapses(dopamine, Excitatory,model='''w:1''',on_pre='''d_post += epsilon_dopa * (dopamine_times(t) * 2 - 1)''',method='exact')
C_I_D = Synapses(Inhibitory, dopamine,model=eqs_S_I, on_pre=eqs_pre_I,method='exact')
C_E_D = Synapses(Excitatory, dopamine,model=eqs_S_E, on_pre=eqs_pre,method='exact')
for A in range(connectivity_matrix.shape[0]):  # Loop through rows
    for B in range(connectivity_matrix.shape[1]):  
        row = connectivity_matrix.index[A]  
        column = connectivity_matrix.columns[B]  
        weight = connectivity_matrix.iloc[A, B]
        weight = float(weight)
        if weight >= 0.0:
            if row == "Excitatory" and column == "Excitatory":
                C_E_E.connect(i=A, j=B)
                C_E_E.w[A, B] = weight
            elif row == "Excitatory" and column == "Inhibitory":
                C_E_I.connect(i=A, j=B)
                C_E_I.w[A, B] = weight
            elif row == "Inhibitory" and column == "Inhibitory":
                C_I_I.connect(i=A, j=B)
                C_I_I.w[A, B] = weight
            elif row == "Inhibitory" and column == "Excitatory":
                C_I_E.connect(i=A,j= B) 
                C_I_E.w[A, B] = weight
            elif row == "Dopamine" and column == "Excitatory":
                C_D_E.connect(i=A,j= B)
                C_D_E.w[A,B] = weight
            elif row == "Dopamine" and column == "Inhibatory":
                C_D_I.connect(i=A,j= B)
                C_D_I.w[A,B] = weight
            elif row == "Inhibatory" and column == "Dopamine":
                C_I_D.connect(i=A,j= B)
                C_I_D.w[A,B] = weight
            elif row == "Excitatory" and column == "Dopamine":
                C_E_D.connect(i=A,j= B)
                C_E_D.w[A,B] = weight
C_E_E.mode = 1
C_E_E.s = '0.01 * rand() * Egmax'
C_E_E.c = 1e-10
C_E_E.d = 0
C_E_I.mode = 1
C_E_I.s = '0.01 * rand() * Egmax'     
C_E_I.c = 1e-10
C_E_I.d = 0
C_I_I.mode = 1     
C_I_I.s = '0.01 * rand() * Igmax'     
C_I_I.c = 1e-10     
C_I_I.d = 0    
C_I_E.s = '0.01 * rand() * Igmax'  
C_I_E.c = 1e-10   
C_I_E.d = 0 
C_I_E.mode = 1
net.add(C_E_E)    
net.add(C_E_I)    
net.add(C_I_I) 
net.add(C_I_E)
synapse_stdp_monitor = StateMonitor(C_E_E, ['s', 'c', 'd'], record=[0])
synapse_stdp_monitor1 = StateMonitor(C_I_E, ['s', 'c', 'd'], record=True, dt=100*ms)
synapse_stdp_monitor2 = StateMonitor(C_E_I, ['s', 'c', 'd'], record=True, dt=100*ms)
synapse_stdp_monitor3 = StateMonitor(C_I_I, ['s', 'c', 'd'], record=True, dt=100*ms)

net.add(synapse_stdp_monitor)
net.add(synapse_stdp_monitor1)
net.add(synapse_stdp_monitor2)
net.add(synapse_stdp_monitor3)

# Stimulate only the inhibitory group
exc_mon = SpikeMonitor(Excitatory)
inh_mon = SpikeMonitor(Inhibitory)
E_mon = StateMonitor(Excitatory, 'v', record= range(103))
v_mon = StateMonitor(Inhibitory,'v', record= range(103))
net.add(exc_mon)
net.add(inh_mon)
net.add(v_mon)
net.add(E_mon)
# Run simulation
net.run(1*second)
# Plot results
figure(figsize=(10, 4))
plot(inh_mon.t/ms, inh_mon.i, 'r.', label="Inhibitory")
xlabel('Time (ms)')
ylabel('Neuron index')
title('Inhibitory Neuron Spikes')

What you have aready tried

Expected output (if relevant)

Actual output (if relevant)

Full traceback of error (if relevant)

2

Hi @asivarajan Could you give some more details on what you were expecting to see and what you are seeing? Unfortunately, I cannot run your code since a number of variables are missing (vr, connectivity_matrix, …).

3

So sorry for the late reply. But I got the code to run and was setting the weights to 0 by accident. The excitatory neurons fire now given that I add an external voltage to a random 10% population. I was wondering on how to add noise as well? Would that be possible?

4

Hi @asivarajan. Glad to hear you solved your issue. There are many ways to add noise in Brian, e.g. you can make your equations stochastic by adding a noise term \xi (see Models and neuron groups — Brian 2 2.7.1 documentation), or you can use PoissonInput to model many incoming spike trains (see Input stimuli — Brian 2 2.7.1 documentation), or you could use functions such as rand or randn in your equations (see Functions — Brian 2 2.7.1 documentation). Let me know if you need help with one of those.

5

Hello @mstimberg thanks so much for your help. I did have another question regarding my connectivity matrix. Right now there should be 47 excitatory neurons, 52 inhibitory neurons, and 4 dopamine neurons. In my connectivity matrix, the total connections are 103 by 103, and the neurons may have matrix indices that are greater the neurongroup indices. For example, there might be an inhibitory neuron in the 78th position in my matrix when theres only indices from 0 to 51. How can i go around that? Right now I have created 103 neurons for each population but there will be neurons that may not be connected.

6

Hi @asivarajan. One thing that helps is to gather all neurons that share the same equations (in your example, excitatory and inhibitory) in a single NeuronGroup, and then use subgroups to distinguish between them:

neurons  = NeuronGroup(num_excitatory + num_inhibitory, eqs, threshold='v>vt', reset='v = vr',
                       method='euler',refractory= 2 * ms)
Excitatory = neurons[:num_excitatory]
Inhibitory = neurons[num_excitatory:]

You can then work with only two instead of four synapse types, since e.g. C_E_E and C_E_I have the same model/effect on the post-synaptic cell:

C_E = Synapses(Excitatory, neurons, model=eqs_S_E, on_pre=eqs_pre, on_post=eqs_post, method='euler')
C_I = Synapses(Inhibitory, neurons, model=eqs_S_I, on_pre=eqs_pre_I, on_post=eqs_post_I, method='euler')

Finally, this allows you to use a column matrix as explained in the documentation (in your example code above you seem to use a DataFrame and not a matrix, but I did not quite understand the exact format; the following assumes a simple weight matrix):

e_sources, e_targets = connectivity_matrix[:num_excitatory, :].nonzero()
C_E.connect(i=e_sources, j=e_targets)
C_E.w[:] = connectivity_matrix[:num_excitatory, :].flatten()
i_sources, i_targets = connectivity_matrix[num_excitatory:num_excitatory+num_inhibitory, :].nonzero()
C_I.connect(i=i_sources, j=i_targets)
C_I.w[:] = connectivity_matrix[num_excitatory:num_excitatory+num_inhibitory, :].flatten()

This is much more efficient than connecting each synapse individually. You can then use the same general approach to create the three types of dopamine connections (dopamine → other neurons, excitatory → dopamine, inhibitory → dopamine). Hope that gets you going!