Hello,
I am trying to model a single LIF neuron. I noticed that I get a different spiking times for Brian2 and my implementation of the model in Matlab (integrating with ODE and in other version solving it “exactly”). An example of spiking times for the same initial conditions is showed below
Brian2 solution:
25.600000000000001, 45.700000000000, 65.79999999999
Exact solution:
25.646636858609416, 45.684117764877215, 65.684799555040613
Maybe someone can explain this behaviour? It seems, that Brian uses time steps equal to 0.1 despite, the “exact” solution method.
The full model description is in https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1002906#s4
Brian Code:
from brian2 import *
import csv
%matplotlib inline
start_scope()
taum = 20*ms
taus= 5 *ms
Vr=-60.
Vth=-40.
freq = 0.05
period = 1./freq*ms
eqs = '''
dv/dt = (Vr-v+I)/taum : 1
dI/dt = -I/taus+mui : 1
mui : Hz
'''
G = NeuronGroup(1, eqs, threshold='v>=Vth', reset='v = Vr', method='exact')
G.v = Vr
G.I = 0.
G.mui = (Vth-Vr)/(taus*(1.-exp(-period/taum)))
M = StateMonitor(G, ['v','I'], record=True)
M2=SpikeMonitor(G,variables='v')
run(500*ms)
savetxt('Vmas.txt',[M.t/ms , M.v[0,:], M.I[0,:]],delimiter=',')
trains = M2.spike_trains();
savetxt('trains.txt',[trains[0]])
Hi. The exact method refers to the fact that we solve the equations exactly in this case instead of using a numerical approximation like forward Euler or Runge-Kutta. This does not change the fact that Brian is a clock-driven simulator, and everything proceeds in discrete time steps. Brian does not support calculating and propagating the exact time of a spike during a time step. The reason is mostly that this is only applicable to a small number of models and that things get either very complicated or very inefficient in large networks where every neuron receives many spikes.
If you have a single neuron (i.e. you don’t care about transmitting the spike somewhere else), there’d be ways to deal with this in Brian by adding a an addition variable hat tracks the time of threshold crossings and an additional term in the reset that corrects for it. But that would be specific to a single model and not applicable in general. Let me know if that interests you and I can give you more pointers.
Thank you for the answer.
My final goal is to model a network of LIF neurons ruled by STDP. As a starting point I was modelling two coupled neurons and comparing results from Brian and Matlab programs. I got different dynamics of synaptic variables (while final states seems to be similar). So my idea was that different dynamics is determined by this discrepancy in spiking times.
I see, yes in the context of STDP the spike timing obviously plays an important role. In case you don’t know the paper: Rudolph&Destexhe discussed this issue in their 2007 article: https://doi.org/10.1016/j.neucom.2006.10.138
There’s one important thing to consider (or potentially to correct if you want to compare things quantitatively): in an exact-timing simulation, you never get the situation that the pre- and the post-synaptic spike occur at exactly the same time. In the clock-driven simulation, however, this can happen, especially with relative big time steps. In the standard formulation of the STDP algorithm with its two trace variables (i.e. as in the Brian example), what effect this has depends on the order in which you handle pre- and post-synaptic spikes. In Brian, the default is to first handle pre-synaptic then post-synaptic spikes, which means this situation (with the default parameters) leads to the maximal possible weight increase. If you evaluated the two in the opposite order, you’d end up with the maximal possible weight decrease. Brian is a bit more explicit about this than other simulators and also allows to change the order (see the scheduling documentation for details). Whatever order is chosen, that also means that the maximal increase and decrease are not exactly symmetric: for one of the two the minimum distance between spikes is 0ms, i.e. it can actually achieve the maximum increase/decrease that was specified; for the other one, the minimum distance between spikes is a single time step, and therefore the maximum increase/decrease is rather 1 - e^{-\frac{\mathrm{dt}}{\tau}} times the specified value. Of course, if either \mathrm{dt} is small or the time constant \tau of your STDP rule is large, this does not make a measurable difference. Still, for exact quantitative comparisons this might be important to keep in mind.
Not sure that was clear at all 
Don’t hesitate to ask if it wasn’t.
Interesting article. I though, that exact timing is not a problem in large networks as “noise” from many inputs or chaotic nature of the model should reduce timing importance. However, in your reference it is showed, that timing affects even spiking rate of all network. I highly doubt that precise timing is very important in real brains, so maybe other formulation of STDP should be used.
I agree that the issue is probably at least partially an artefact of the way we describe STDP. Numerical accuracy and related issues are certainly important (e.g. if your simulation gives qualitatively different results with dt=0.1ms and dt=0.05ms this is worrying), but I feel that we sometimes focus too hard on precise simulation of crude models…