Izhikevich neuron model : is my understanding of the reset of the variable "u" wrong?

Support
1

Description of problem

Hi,

In the Izhikevich model, when a spike occurs, i.e the membrane potential reached the threshold voltage, its variables “v” (membrane potential) and “u” (membrane recovery) are reset as :

v = c
u = u + d

In BRIAN, I noticed that when a spike occurs, the variables reset (in discrete) looks like this :

v[n+1] = c
u[n+1] = u[n] + d + du

And so, it’s actually the “du” expression that disturbs me. Is it not supposed to be just like this ?

v[n+1] = c
u[n+1] = u[n] + d

Depending on the timestep value, the emission of spikes with constant current stimulation is different when comparing both resetting methods.

Minimal code to reproduce problem

This is my Brian code that simulates one IZH neuron and displays its variables “v” and “u” at each timestep until the first spike.

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from brian2 import *

start_scope()

duration = 0.03*second          # simulation time to observe the first spike 
defaultclock.dt = 1*ms           # time step dt 

# IZH parameters
I_exc = 7.0
a = 0.02
b = 0.2 
c = -65.0
d = 8.0
v_th = 30.0

v_init = c
u_init = b*c

time_refractory = 0.0*ms

# IZH equations 
eqs_IZH = '''
        	dv/dt = (0.04*v**2 + 5*v + 140 - u + I_exc)/ms : 1
       	 	du/dt = a*(b*v - u)/ms : 1
      	  '''
      	 
reset = '''
          v = c
          u += d
        '''    

# 
NG = NeuronGroup(1, eqs_IZH, threshold='v>v_th', reset=reset, refractory=time_refractory, method="euler")
NG.v = v_init
NG.u = u_init

# MONITORS
s_mon = SpikeMonitor(NG)
v_mon = StateMonitor(NG, variables = True, record = True)


# RUN SIMULATION
run(duration)



# DEBUG : printing variables "v", "u" and emitted "spike" for each timestep  
time_spike = 0
print ('////////////////////////////////////////////////////')
print("t \t\t  v \t\t   u \t\t spike")
for i in range(0,int(s_mon.t[0]/defaultclock.dt)+3):   # print until the 1st spike + 2 timesteps after
	if (i == int(s_mon.t[time_spike]/defaultclock.dt)):	
		print("%2.6f \t %2.6f \t %2.6f \t 1 <-----" % (v_mon.t[i], v_mon.v[0][i], v_mon.u[0][i]))
		if (time_spike != size(s_mon.t)-1):
			time_spike = time_spike + 1
	else :
		print("%2.6f \t %2.6f \t %2.6f \t 0" % (v_mon.t[i], v_mon.v[0][i], v_mon.u[0][i]))
print ('///////////////////////////////////////////////////')

Expected output (if the proposed “u” reset expression is correct)

Here I manually added the expected values of “u” after a spike. Therefore, the next spike occurs at 0.061 ms for dt=1ms.

////////////////////////////////////////////////////
t            v 	              u           spike
0.000000 	 -61.000000 	 -13.000000 	 0
0.001000 	 -57.160000 	 -12.984000 	 0
0.002000 	 -52.285376 	 -12.952960 	 0
0.003000 	 -44.408874 	 -12.903042 	 0
0.004000 	 -27.664279 	 -12.822617 	 0
0.005000 	 24.449437 	     -12.676822 	 0 #    u[n+1]       u[n]     d      
0.006000 	 -65.000000 	 -4.676822	     1 #<----- u   = -12.676822 + 8 
0.007000 	 -69.3232... 	 -4.8432...	     0
0.008000 	 -71.8677... 	 -5.0237...	     0
0.009000 	 -72.5838... 	 -5.2107...      0
[...]           [...]           [...]     [...] # manually added
0.061000     -65.0000000     -1.9779...      1  # next spike
///////////////////////////////////////////////////

Actual output

As above, for readability, I manually added the time of the next spike at 0.062 ms for dt=1ms.

////////////////////////////////////////////////////
t            v 	              u           spike
0.000000 	 -61.000000 	 -13.000000 	 0
0.001000 	 -57.160000 	 -12.984000 	 0
0.002000 	 -52.285376 	 -12.952960 	 0
0.003000 	 -44.408874 	 -12.903042 	 0
0.004000 	 -27.664279 	 -12.822617 	 0
0.005000 	 24.449437 	     -12.676822 	 0  #    u[n+1]       u[n]     d      du
0.006000 	 -65.000000 	 -4.325488 	     1  #<----- u   = -12.676822 + 8 + 0.351334...
0.007000 	 -69.674512 	 -4.498978 	     0
0.008000 	 -72.366589 	 -4.687696 	     0
0.009000 	 -73.034910 	 -4.883409 	     0
[...]           [...]           [...]     [...] # manually added
0.062000     -65.000000      -1.922917       1  # next spike
///////////////////////////////////////////////////

So is my understanding of the reset of the variable “u” wrong ?

Thanks for your help,

Marino

2

Hi Marino. Just to make sure I understand correctly, with du you mean the update term (e.g. dt\cdot f'(u) when using the Forward Euler algorithm)? The mathematical formulation of the Izhikevich model does not specify this exactly, but I think the standard Brian implementation is the usual one. If you wanted to update u based only on the value of the previous time step, then you’d have to split the update of v and u, i.e. do something like this:

new_v = v[n] + dv
if new_v > v_threshold:
   v[n+1] = c
   u[n+1] = u[n] + d
else:
   v[n+1] = new_v
   u[n+1] = u[n] + du

The standard implementation rather updates all state variables together:

new_v = v[n] + dv
new_u = u[n] + du  # note that du depends on v[n], so we can't overwrite it yet
v[n+1] = new_v
u[n+1] = new_u
if v[n+1] > v_threshold:
  v[n+1] = v
  u[n+1] = u[n+1] + d

In principle, you could implement the first approach in Brian as well, but it would be a bit cumbersome. Also, Izhikevich’s own code uses approach 2 (see e.g. https://www.izhikevich.org/publications/net.m), but with the state update after threshold check and reset (i.e. the spike is emitted in the time step after the threshold crossing).

1 Like
3

Hi Mstimberg,

Yes you are correct, I call du the update term of the u variable when using the Forward Euler algorithm.

I’m taking BRIAN’s simulations as a reference to my hardware implementations. And now, I realised that my hardware resetting implementation of the izhikevich model (the first you described) is not correct.

There are different implementations outside there :confused: But mostly are the same as BRIAN (so following the Izhikevich paper way), and it makes sense to me now.

So thank you again for your time and making it clearer !