Firing rate of individual neurons

Support
1

Description of problem

Hello Brian2 team,

I am just wondering if you have any insight into what is the best way to get an array of firing rates of individual neurons within a network. So if I have a population of 50 neurons and I simulated it for 1000 time steps, I’d want an array of rates of size 50 x 1000?

What you have already tried

I have used the PopulationRateMonitor but it seems silly to create a monitor for each neuron? Does it make more sense to use the SpikeMonitor and process the output times and indices manually to create the individual rates?

Sorry if I’m missing something obvious! Thanks for your help =)

2

HI Imun,

My usual solution for this is to create a single spike monitor for the whole population, then bin each neuron’s spike times to get a vector of n time_steps with an integer for spike count, then apply a gaussian filter to smooth this vector into an estimate of rate.

The other solution I’ve tried which is less flexible, but gives you an online estimate of firing rate, and can be monitored simply with a StateMonitor is to set up an exponential filter “neuron” connected to each neuron you want to record. The downsides are this estimate can only look at preceding spikes so it will generally be a less smooth estimate, and also your choice of time-constant for filtering the rate has to be decided before the simulation, so it’s not easy to adapt the smoothing based on your observations after the simulation.

It sounds like you’ve already considered the first of these approaches, but let me know if you have trouble implementing either and I can help.

Otherwise, I think there’s a sensible feature request in there:

1 Like
3

Thanks Adam! This was helpful =) I managed to create the frequency matrix I needed for the jPCA function to show low dimensional rotational dynamics in my neuron groups. Here is my code (might be quite naive so feel free to let me know if there are more efficient ways) and a cool image to share with the community.

The dynamical solution also sounds interesting and potentially useful for me in the future so I’ll keep this in mind! The feature request would also be cool for the future. Thanks again.

  num_neurons = 500 
        time_vector=np.arange(0, 300, 5)
        frequency_matrix = np.zeros([len(time_vector), num_neurons])
        times = spikemon.t/b2.ms
        indices = spikemon.i
        timestep = 0.005 #ms
        for i in range(len(time_vector)-1):
            time_mask = ((times >= time_vector[i]) & (times <= time_vector[i + 1]))
            indexbin = indices[time_mask]
            for j in range(num_neurons):
                num_count=indexbin.tolist().count(j)
                frequency_matrix[i][j] = num_count/timestep
                
        for k in range(frequency_matrix.shape[1]):
            frequency_matrix[:, k] = scipy.ndimage.gaussian_filter1d(frequency_matrix[:, k], sigma=6)

qhFPoqLrpV

1 Like
4

Thanks for sharing the code (and the cool image :slight_smile: )! There’s some numpy magic that you can use to calculate the frequency matrix in a very fast way (though your current solution is probably fast enough for most purposes):

frequency_matrix = np.zeros([len(time_vector), num_neurons])
# Add 1 in the correct place for each spike:
np.add.at(frequency_matrix,
          (np.array(spikemon.t / (5*b2.ms), dtype=int), spikemon.i), 1)
frequency_matrix /= timestep

I think this gives the exact same values as your solution, except for some numerical issues: A spike that occurs exactly at, say, 10*ms might end up in the 5ms–10ms or in the 10ms–15ms bin.

I agree! This would probably have to be in SpikeMonitor and not in PopulationRateMonitor, though, since the latter does not store the information for the individual neurons.

5

just submitted as a feature request on github:

1 Like