LIF spiking times

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.