Optimization of simulation time using TimedArray

Support
1

Description of problem

Dear all, I’m trying to simulate a neuron population receiving an external input. Specifically, neurons are all-to-all connected and have one state variable phi evolving according to a differential model. In addition, I have an external input which is a 1-dimensional variable A with a known value at every time instant.

To distribute this input to all neurons, at every integration step t I want to sample one value per neuron from a normal distribution centered on A(t). One sampled value enters each neuron’s differential equation additively (see eq. below).

To implement this, I first create a TimedArray of size NxT (n. neurons, n. time steps), which stores all sampled values to feed neurons at every integration step. The resulting simulation is (to the best of my knowledge) quite slow, does anybody have a better implementation to propose?

Minimal code to reproduce problem

import brian2 as br
dtime = 100*br.ms
T = 2000*br.ms
N = 5000
A = br.random.rand(int(T/dtime)) # test input

# configure
br.start_scope()
br.defaultclock.dt = dtime

# define external drive
# sample inputs to each neuron from a normal distribution centered in A[t]
samples = br.array([br.random.normal(size=N) for i in range(A.shape[0])])
A = br.TimedArray(br.tile(A,(samples.shape[1],1)).T * samples,dt=dtime)

# define equations
eqs = '''
dphi/dt = omega_int + a_int*sin(phi) + coupl + cos(phi)*A(t,i)*Hz + sigma_int*xi : 1
coupl : hertz
omega_int : hertz
a_int : hertz
sigma_int : hertz**0.5
'''

# define population
in_cond = 2*math.pi * br.np.random.uniform(size=(N))
neurons = br.NeuronGroup(N,eqs,method='heun')
neurons.phi = in_cond
neurons.omega_int = 1*br.Hz
neurons.a_int = 1*br.Hz
neurons.sigma_int = 1*br.Hz**0.5

# define connections
J = 1*br.Hz
conn = br.Synapses(neurons,neurons,'coupl_post = J/N_pre*sin(phi_pre-phi_post) : hertz (summed)')
conn.connect(p=1)

# define monitor
activity = br.StateMonitor(neurons,'phi',record=True,dt=dtime)

# run
br.run(T)

What you have aready tried

I have some ideas on what might be slowing down my simulation:

  • Creating the whole TimedArray before simulating requires a slow access to each of its elements during the simulation?
  • A StateMonitor recording a variable at every timestep is a slow operation?

I think I need a more educated opinion on how to efficiently use brian2.

Expected output (if relevant)

Simulating the same dynamical system, optimizing simulation time.

Actual output (if relevant)

Full traceback of error (if relevant)

2

Hi @pietro. In this case, the TimedArray does not seem to affect the simulation time – if you replace the A(t, i) term in your equations by a fixed value, the simulation time will not change. The simulation runs slow because of the summed variable operation that updates the coupling term every timestep for all 25 000 000 connections. If you’d be ok with updating this term only every second time step for example, the time would drastically go down. Brian has a built-in profiling mechanism that you can use to figure out where time is spent. In this case, it is almost 100% in the summed variable update:

# run
br.run(T, report='text', profile=True)
print(br.profiling_summary())

gives

Starting simulation at t=0. s for a duration of 2. s
2. s (100%) simulated in 6s
Profiling summary
=================
synapses_summed_variable_coupl_post    6.63 s    99.92 %
neurongroup_stateupdater               0.00 s     0.06 %
statemonitor                           0.00 s     0.02 %
neurongroup                            0.00 s     0.00 %
synapses                               0.00 s     0.00 %

Hope that clears things up!