I would like to create a SpikeGeneratorGroup with multiple “neurons” in it each running at a different frequency. I can specify the spike times for each neuron given their frequency but given the length of simulation and number of frequencies it on the order of 0.5 billion time points which is not ideal.
Alternatively I’ve tried creating a number of SpikeGeneratorGroup objects each with one neuron and connecting them independently and after about ~20 or so synapses it takes longer and longer to create the synapses (bordering on not feasible for larger numbers).
Is there a better/good way to do this?
Code to generate spike times
N = 20
neurons = NeuronGroup(N, model='model code here')
T = [10, 20, 30, 40, 50] # period in ms
times = [np.arange(0, 1000, per) for per in T]
ixs = [[ix] * len(s) for ix, s in enumerate(times)]
times = np.concatenate(times) * ms
ixs = np.concatenate(ixs)
input = SpikeGeneratorGroup(N, ixs, times) # periodic input @ 400Hz
synapses = Synapses(input, neurons, on_pre='update code here') # connect input to neurons
synapses.connect('i==j') # one synapse goes to one neuron
Code to generate multiple objects
N = 20
neurons = NeuronGroup(N, model='model code here')
input = []
synapses = []
for ix in range(N):
input.append(SpikeGeneratorGroup(1, [0], [0] * ms, period= N * ms))
synapses.append(Synapses(input[-1], neurons, on_pre="update code here"))
synapses[-1].connect(i=0, j=ix)
1 Like
Great question! I agree that it does not make much sense to precalculate all the spike times, given how regular they are. In your second approach, I suspect it is not really the synapse creation that takes a long time, but rather the generation/compilation of the code. In general, having many separate simple objects is very inefficient in Brian.
A better solution would be to abandon SpikeGeneratorGroup and instead use a NeuronGroup that calculates the spike times online.
One approach would basically model them as a crude (non-leaky) integrate-and-fire neuron with a constant input current that triggers regular spiking with the intended rate:
input = NeuronGroup(N, '''dx/dt = rate : 1
rate : Hz (constant)''',
threshold='x>1', reset='x=0')
A slightly more abstract/efficient solution would directly check times:
input = NeuronGroup(N, '''period : second (constant)
lastspike : second''',
threshold='(t-lastspike)>period', reset='lastspike=t')
Note that in both cases, spike times might shift a bit during very long runs due to rounding issues. If this is a problem you could instead pre-calculate an integer period (in time steps), and then trigger spikes based on the t_in_timesteps variable (e.g. with an integer_period : integer (constant) variable replacing period, and then checking t_in_timesteps % integer_period == 0 as the threshold condition).
Thank you for that quick reply, That’s a really elegant solution, I quite like it.
If you were to guess, how long would long be before drift starts to become an issue? We’re only looking in the ~2 second simulation window.
I just realized that this is actually not about drift (because we are resetting things regularly) but rather about the imprecision when comparing things. If the period is a multiple of Brian’s simulation time step, then we are supposed to land exactly on the threshold value which means in practice we are either very slightly below or above it which will change the period by a single time step. So for best results I’d either use the integer approach, or the time-based solution but with Brian’s robust conversion to integer timesteps:
input = NeuronGroup(N, '''period : second (constant)
lastspike : second''',
threshold='timestep(t-lastspike, dt)==timestep(period, dt)', reset='lastspike=t')
But maybe it is worth running only the input first, recording it with a spike monitor and checking whether the diff of the spike times is what you expect.