A way to save the brain2 objects for later use + single neuron average FR

Support
1

Description of problem

I have a simulation which takes a while to run. Even after parallelisation. Now I was wondering if there is a way to save the array exactly as to be saved in data and import and use later. What is the best way to do this ?

Plus :

Is there an efficient way to average the single neuron firing rate over n_sims. that is , neuron1 average the firing rate 1000 times . I could’nt setup the PopulationMonitor rate to do so. I realise maybe its is not the thing to use for single neurons.

Minimal code to reproduce problem

#imports
from brian2 import *
from IPython.display import clear_output
import numpy as np
import matplotlib.pyplot as plt
set_device('cpp_standalone', build_on_run=False)
prefs.devices.cpp_standalone.openmp_threads = 12

#define ovrlap AUB etc
n_ab = 0.33 # this percentage receives both sine and pulse
n_a = (1 - n_ab)/2 # only sine
n_b = (1 - n_ab)/2 # only pulse


########## neuron inputs #############
start_scope()
defaultclock.dt = 0.05*ms
duration = 2*second

Cm = 1*uF  # /cm**2
gL = 0.1*msiemens
EL = -65*mV
ENa = 55*mV
EK = -90*mV
gNa = 35*msiemens
gK = 9*msiemens

# Synaptic variables
taue = 5*ms
taui = 10*ms
Ee = 0*mV
Ei = -80*mV
we = 15*uS  # excitatory synaptic weight
wi = 67*uS  # inhibitory synaptic weight

# Hodgkin-Huxley + synaptic + external currents
eqs = '''
dv/dt = (-gNa*m**3*h*(v-ENa) - gK*n**4*(v-EK) - gL*(v-EL)
         + Iamp0 + Iamp * sin(2*pi*Ifreq*t)
         + ge*(Ee-v) + gi*(Ei-v)) / Cm : volt

# gating vars
m = alpha_m/(alpha_m+beta_m) : 1
alpha_m = 0.1/mV*10*mV/exprel(-(v+35*mV)/(10*mV))/ms : Hz
beta_m = 4*exp(-(v+60*mV)/(18*mV))/ms : Hz

dh/dt = 5*(alpha_h*(1-h)-beta_h*h) : 1
alpha_h = 0.07*exp(-(v+58*mV)/(20*mV))/ms : Hz
beta_h = 1./(exp(-0.1/mV*(v+28*mV))+1)/ms : Hz

dn/dt = 5*(alpha_n*(1-n)-beta_n*n) : 1
alpha_n = 0.01/mV*10*mV/exprel(-(v+34*mV)/(10*mV))/ms : Hz
beta_n = 0.125*exp(-(v+44*mV)/(80*mV))/ms : Hz

# synaptic conductances
dge/dt = -ge*(1./taue) : siemens
dgi/dt = -gi*(1./taui) : siemens

# external input parameters
Iamp0 : amp   # constant or pulsed current
Iamp : amp    # sinusoidal amplitude
Ifreq : Hz    # sinusoidal frequency
'''

# Neuron group
N = 100
group = NeuronGroup(N, eqs, threshold='v>-20*mV', refractory=3*ms,
                    method='rk4')
group.v = -70*mV
group.h = 1


######### noise goes everywhere ##########
# Noise input
rate = 20*Hz
Noise_in = PoissonGroup(N, rates=rate)
S = Synapses(Noise_in, group, on_pre='ge += we')
S.connect(j='i')

# Example: set input currents in n_a neurons

n_a_arr = np.zeros(N)
# first n_a percentage becomes 1
n_a = n_a * N  # convert percentage to number of neurons
n_a_arr[:int(n_a)] = 1  # set first n_a neurons to 1  

# sinusoidal input
n_b_arr = np.zeros(N)
#last n_b percentage becomes 1
n_b = n_b * N  # convert percentage to number of neurons
n_b_arr[-int(n_b):] = 1  # set last n_b neurons to 1

n_ab_arr = n_a_arr + n_b_arr  # overlap neurons
#make 0s as 1s and 1s as 0s
n_ab_arr[n_ab_arr == 0] = 1
n_ab_arr[n_ab_arr == 1] = 0

# Monitors
monitor = StateMonitor(group, 'v', record=True)
monitor_s = SpikeMonitor(group)
# syn_rate = PopulationRateMonitor(Noise_in)
# monitor_rate = PopulationRateMonitor(group)
# ge_monitor = StateMonitor(group, 'ge', record=True)
# input_monitor = SpikeMonitor(Noise_in)


n_sims = 1000

from tqdm import tqdm
for i in tqdm(range(n_sims)):    
    group.Iamp0 = 0*0.55*uA * (n_a_arr + n_ab_arr)   # baseline or pulse amplitude
    group.Iamp  = 0*0.35*uA * (n_b_arr + n_ab_arr)  # sinusoidal amplitude
    group.Ifreq = 8*Hz     # sinusoidal frequency

    run(duration/5.0)

    group.Iamp0 = 0.55*uA * (n_a_arr + n_ab_arr)   # baseline or pulse amplitude
    group.Iamp  = 0.35*uA * (n_b_arr + n_ab_arr)  # sinusoidal amplitude
    group.Ifreq = 8*Hz     # sinusoidal frequency

    run(3*duration/5.0)

    group.Iamp0 = 0*0.55*uA * (n_a_arr + n_ab_arr)   # baseline or pulse amplitude
    group.Iamp  = 0*0.35*uA * (n_b_arr + n_ab_arr)  # sinusoidal amplitude
    group.Ifreq = 8*Hz     # sinusoidal frequency   

    run(duration/5.0)

device.build()
clear_output(wait=False)
print('Whole sim completed')

#plot the raster plot for neuron group
plt.figure(figsize=(10, 6), dpi=150)
plt.plot(monitor_s.t/ms, monitor_s.i, 'k|', markersize= 3)
plt.xlabel('Time (ms)', fontsize=12)
plt.ylabel('Neuron index', fontsize=12)
plt.title('Raster plot of neuron spikes', fontsize=14)
plt.tight_layout()
plt.show()


#plot the raster plot for neuron group
plt.figure(figsize=(10, 6), dpi=150)
plt.plot(monitor_s.t/ms, monitor_s.i, 'k|', markersize= 3)
plt.xlabel('Time (ms)', fontsize=12)
plt.ylabel('Neuron index', fontsize=12)
plt.title('Raster plot of neuron spikes', fontsize=14)
plt.tight_layout()
plt.xlim(0, 4000)
plt.show()

What you have already tried

saving the brian2 object : have’nt tried anything yet.

For single neuron type firing rate average : I have tried chopping off the numpy array with time stamps. For example, unpacking the whole spike time and subtracting 2000, 4000, 6000 etc and pool them. But it is very messy.

Expected output (if relevant)

an array of Nxm , where N = 100 (number of neuron) x m = int(2000/dt) (number of time steps). Each row giving the Average firing rate for single neuron

Actual output (if relevant)

Full traceback of error (if relevant)

2

Hi @Jeevan_codes8580

That’s a great question. Running many trials and analyzing the results efficiently is a common challenge … I am more familiar with the code side of brain2 and don’t have much hands on experience writing and running models … but what I can suggest is this

1. Saving and Loading Simulation Data

For cpp_standalone mode, the most robust and efficient way to save your data is to collect it during the simulation and then save the raw NumPy arrays. The network.store() and network.restore() methods are not implemented for this mode ( if I am correct) , but saving the data directly with NumPy is straightforward and effective.

Here’s how you can structure your script to save the results after the simulation is complete:

# After your simulation loop finishes...

# 1. Collect all the data you want to save
# For example, the total spike counts and the raster plot of the final trial
simulation_data = {
    'total_spike_counts': total_spike_counts,
    'final_trial_t': final_trial_t,
    'final_trial_i': final_trial_i,
}

# 2. Save the data to a single, compressed file
# This is efficient and keeps all your results together
np.savez_compressed('simulation_results.npz', **simulation_data)

print("Simulation data saved to simulation_results.npz")

# 3. To load the data back in a separate script or later on:
loaded_data = np.load('simulation_results.npz')

# You can now access the arrays like a dictionary
# loaded_data['total_spike_counts'], loaded_data['final_trial_t'], etc.

This method is fast, creates a single organized file …

2. Calculating Per-Neuron Firing Rates

You are right, PopulationRateMonitor gives the average for the whole group. To get the average rate for each individual neuron across all your trials, the best tool is the SpikeMonitor.

what you can do is something like :

import numpy as np

def calculate_firing_rates_per_neuron(spike_monitor, n_neurons, duration, n_sims):
    """
    Calculates the average firing rate for each neuron across multiple trials.

    Parameters:
    - spike_monitor: The Brian2 SpikeMonitor from the final trial.
    - total_spike_counts: A numpy array with the total spike count for each neuron.
    - n_neurons: The total number of neurons (N).
    - duration: The duration of a single simulation trial (as a Brian2 Quantity).
    - n_sims: The total number of simulations run.

    Returns:
    - average_firing_rates: A Brian2 Quantity array with the average firing rate in Hz.
    """
    # This array should be accumulated over all trials
    total_spike_counts = np.bincount(spike_monitor.i, minlength=n_neurons)

    # Calculate average rate: total spikes / total time
    total_duration = n_sims * duration
    average_firing_rates = total_spike_counts / total_duration

    return average_firing_rates

# --- USAGE IN YOUR SCRIPT ---
# This would go inside your simulation loop to accumulate counts,
# and the final calculation would be after the loop.

# (Like Inside the loop, after each run, do this )
# spike_counts_this_trial = np.bincount(monitor_s.i, minlength=N)
# total_spike_counts += spike_counts_this_trial

# (and then after the loop)
# average_rates = total_spike_counts / (n_sims * duration)
# print(f"Average firing rates: {average_rates}")
# print(f"Mean across all neurons: {np.mean(average_rates):.2f}")

This approach is much faster than iterating through each neuron, especially for a large number of neurons or spikes.

By the way, keep an eye on developments for RateMonitor we are working on here Feat: Introduce abstract RateMonitor class for unified rate analysis by Legend101Zz · Pull Request #1657 · brian-team/brian2 · GitHub , once that’s merged , the new binned() method will make calculating time-binned population rates more flexible in the future, but for per-neuron rates across trials, the SpikeMonitor method is the way to go for now , I guess …

1 Like
3

thanks this is quite helpful. I will keep an eye on ratemonitor

1 Like
4

Hi all. Just to add to this discussion: two things that might be helpful (I have to admit, I did not 100% understand what you want your average spike rates to be averaged over – trials? or over time windows within a single simulation?):

  • The SpikeMonitor already has a count attribute, which is an array of the total spike count for each neuron, i.e. if you divive this array by the simulation time, you have an array of the average firing rate of each neuron
  • The SpikeMonitor stores time/indices for each spike, but you can convert this into a matrix of 0/1 values for each neuron and time step, which might make later processing easier:
    indices, times = monitor_s.it
matrix = np.zeros((len(group), int(round(duration/defaultclock.dt))), dtype=bool)
matrix[times, indices] = True
    For example, you can then get the spike count for non-overlapping windows of W timesteps with reshaping and summing (note that this only works if the total number of timesteps is a multiple of W):
    binned = matrix.reshape((len(group), -1, W)).sum(axis=-1)