Two-Layer Network in Brian2: Issue with Firing Rate and Spiking Activity

Hey everyone,

I’m working on a two-layer network in Brian2 and encountering some strange behavior. The first layer of neurons receives a fixed input rate for each neuron, and I’m defining the NeuronGroup like this:

globals()[name_prefix] = NeuronGroup(total_neurons, inp_eqs, threshold='rand() < rate*dt',
                                    method='euler', name=name_prefix)

The input equation (inp_eqs) is as follows:

inp_eqs = '''
    unit_idx : integer (constant)
    x : 1
    y : 1
    rate = input_array(t, i) : Hz
'''

For the input, I want the neurons to receive a specific rate for the first 100ms, and then for the next 100ms, only the background firing should be present. The input array is constructed as follows:

# Combine input from two populations with their respective background rates
input_array = np.hstack((self.amp_inp[0] * input[0] + self.bg_inp_rate[0],
                         self.amp_inp[1] * input[1] + self.bg_inp_rate[1]))

# Create an array of background rates for the second 100ms
half_length = input_array.size // 2
zero_array = np.hstack((np.full(half_length, self.bg_inp_rate[0]),
                        np.full(half_length, self.bg_inp_rate[1])))

# Stack the input rates and background rates
combined_array = np.vstack([input_array, zero_array])

# Use TimedArray to apply this input over time
input_mat = np.tile(combined_array, (1, 1))  # Adjust the tiling if needed
self.input_array = TimedArray(input_mat * Hz, dt=self.dur_rep)

After this first layer, the second layer is receiving the spiking activity of the neurons in the excitatory population and connects to an inhibitory population. However, I’m seeing some weird spiking activity in the second layer.

Problems I’m facing:

  1. The firing rate of the first layer neurons doesn’t seem to behave as expected. I’m seeing odd fluctuations that don’t match the input rate over time.
  2. The second layer, which should be influenced by the first layer’s spiking activity, shows strange, unpredictable patterns of spiking, especially high spiking activity at the begining of the simulation.

Does anyone have any ideas about why this might be happening? Any tips on debugging the firing rates or dealing with spiking activity transfer between layers in Brian2 would be much appreciated.

Thanks in advance!


Hi @elnaz91. I don’t see anything that is obviously wrong, but the devil is often in the details. Regarding the TimedArray, maybe verify that the shape of the array is correct? It should be timesteps × neurons. I don’t completely follow the creation of input_mat in your code – when you talk about “two populations”, do you mean that you only have two different inputs? In the equations, you use input_array(t, i), so you need to provide one input per neuron (or use something like input_array(t, i%2) or input_array(t, i//(N/2)) to provide two different input types to the neurons). That said, if your input has only two possible values for each neuron (background or background + fixed rate), then it would probably easier to define this directly in the equations. But not sure I understood correctly.

Regarding the strong activity at the beginning: this is often due to a missing initialisation. If your neuron model for the second layer uses “biological values” (e.g. say -50*mV for the threshold and not an abstract threshold such as 1), then the initial value should not be 0 (the default), because otherwise all neuron spike in the initial time step. If you already initialized the membrane potential (and potentially other values in a complex model), then the problem might be that they at the same level for all neurons, consider using a random initialisation instead.

Hope that gives you a few ideas, if this doesn’t work out, do not hesitate to share more details!

1 Like

Hi Marcel, Thank you so much for your help! The issue with strong activity has been resolved by adjusting v_int.

Now, for the first problem, let me clarify what I’m trying to achieve:

I’m working with image patches of size 16x16. After flattening and processing them, I convert these patches into on/off states. These on/off patches represent the rates fed into my network. Since the size is 256, I need 512 neurons to maintain a 1:1 correspondence.

My goal is to implement an interval representation of the input. I present the patches to the network for 100 ms, followed by a rest period,. However, I am a bit confused about how to use the Timed Array feature—I’m not sure I fully understand how it works.

Additionally, I would prefer to define neuron groups for the network rather than using Poisson neurons.

P.S The the code that I shared previously is how I initiazlized my input…

        #####################################################################################
        ##                                  Layer Input                                    ##
        #####################################################################################
        if ly_type == 'input':
            unpack_class_attributes(self, globals())
            total_neurons = 2*self.Nunits * self.Nneuron_input

            globals()[name_prefix] = NeuronGroup(total_neurons, inp_eqs, threshold='rand() < rate*dt', 
                               method='euler', name=name_prefix)
            globals()[name_prefix].x        = '(i % (Munits*Mneuron_input)) * (1.0/Mneuron_input)'
            globals()[name_prefix].y        = '(i // (Munits*Mneuron_input)) * (1.0/(2*Mneuron_input))'
            globals()[name_prefix].unit_idx = 'floor(x) + floor(y)*Munits'
                      
    input_params = {   
    'inp_eqs': '''
        unit_idx : integer (constant)
        x : 1
        y : 1
        rate = input_array(t, i) : Hz 
    ''',    
    'dur_rep'     : 50*ms,
    'bg_inp_rate' : [17, 8],
    'amp_inp'     : [1e+2, 1e+2]
}
    hsnn.set_parameters(input_params)
    hsnn.network_input(input)
    def network_input(self, input):
        input_array = np.hstack((self.amp_inp[0] * input[0] + self.bg_inp_rate[0], self.amp_inp[1] * input[1] + self.bg_inp_rate[1]))
        half_length = input_array.size // 2  # Assuming input_array has an even number of elements
        zero_array = np.hstack((np.full(half_length, self.bg_inp_rate[0]),
                                np.full(half_length, self.bg_inp_rate[1])))
        combined_array = np.vstack([input_array, zero_array])
        input_mat = np.tile(combined_array, (2, 1))  # Atwo times repeats

        self.input_array = TimedArray(input_mat * Hz,  dt = self.dur_rep)

This looks reasonable to me, but what is the shape of your input array. From what you say, it should have shape (2, 256), i.e. the on and off values for each pixel. Your code seems to use

50ms stimulus - 50 ms silence - 50ms stimulus - 50ms silence
(with the same stimulus), but I don’t think this is what you are worrying about?

I used a variant of your code to plot the input for three neurons, and it looks reasonable, I’d say. Or did you expect something else?

from brian2 import *
import numpy as np
amp_inp = [100, 100]
bg_inp_rate = [17, 8]
dur_rep = 50*ms
input = np.random.rand(2, 256)

input_array = np.hstack((amp_inp[0] * input[0] + bg_inp_rate[0], amp_inp[1] * input[1] + bg_inp_rate[1]))
half_length = input_array.size // 2  # Assuming input_array has an even number of elements
zero_array = np.hstack((np.full(half_length, bg_inp_rate[0]),
                        np.full(half_length, bg_inp_rate[1])))
combined_array = np.vstack([input_array, zero_array])
input_mat = np.tile(combined_array, (2, 1))  # Atwo times repeats
input_array = TimedArray(input_mat * Hz,  dt = dur_rep)

times = np.arange(200)*ms
# plot input for a few on neurons
for idx in range(3):
    plt.plot(input_array(times, idx))
for idx in range(3):
    plt.plot(input_array(times, 256 + idx), ls='--')
plt.show()

But the way you use the TimedArray might be a bit too complex here. You could instead use the TimedArray only to denote your stimulus/silence times and assign the input directly to your group. I.e. something along the lines of

input_eqs = '''
unit_idx : integer (constant)
        x : 1
        y : 1
        input : Hz
        background : Hz
        rate = background + input * stimulus_on(t) : Hz 
'''
stimulus_on = TimedArray([1, 0, 1, 0], dt=50*ms)
input_group = NeuronGroup(512, ...)
on_group = input_group[:256]
off_group = input_group[256:]
on_group.background = bg_inp_rate[0]
on_group.input = input[0]
off_group.background = bg_inp_rate[1]
off_group.input = input[1]

Hope that helps!