COMMENT
A "simple" implementation of the Izhikevich neuron.
Equations and parameter values are taken from
Izhikevich EM (2007).
"Dynamical systems in neuroscience"
MIT Press
Equation for synaptic inputs taken from
Izhikevich EM, Edelman GM (2008).
"Large-scale model of mammalian thalamocortical systems."
PNAS 105(9) 3593-3598.
Example usage (in Python):
from neuron import h
sec = h.Section(name=sec) # section will be used to calculate v
izh = h.Izhi2007b(0.5)
def initiz () : sec.v=-60
fih=h.FInitializeHandler(initz)
izh.Iin = 70 # current clamp
Cell types available are based on Izhikevich, 2007 book:
1. RS - Layer 5 regular spiking pyramidal cell (fig 8.12 from 2007 book)
2. IB - Layer 5 intrinsically bursting cell (fig 8.19 from 2007 book)
3. CH - Cat primary visual cortex chattering cell (fig 8.23 from 2007 book)
4. LTS - Rat barrel cortex Low-threshold spiking interneuron (fig 8.25 from 2007 book)
5. FS - Rat visual cortex layer 5 fast-spiking interneuron (fig 8.27 from 2007 book)
6. TC - Cat dorsal LGN thalamocortical (TC) cell (fig 8.31 from 2007 book)
7. RTN - Rat reticular thalamic nucleus (RTN) cell (fig 8.32 from 2007 book)
ENDCOMMENT
: Declare name of object and variables
NEURON {
POINT_PROCESS Izhi2007b_dyn_thr
RANGE C, k, vr, v0, u0, vt, vt0, vpeak, dyn_th, u, a, b, c, d, alpha, Vi, Vmin, Ka, Ki, tau, Iin, celltype, alive, cellid, verbose, derivtype, delta, t0
NONSPECIFIC_CURRENT i
}
: Specify units that have physiological interpretations (NB: ms is already declared)
UNITS {
(mV) = (millivolt)
(uM) = (micrometer)
}
: Parameters from Izhikevich 2007, MIT Press for regular spiking pyramidal cell
PARAMETER {
C = 1 : Capacitance
k = 0.7
vr = -60 (mV) : Resting membrane potential
v0 = -70 (mV) : Membrane potential (initial condition)
u0 = 0 : Recovery variable (initial condition)
vt0 = -40 (mV) : Membrane threhsold (initial condition, if needed)
vpeak = 35 (mV) : Peak voltage
dyn_th = 0 : dynamic threshold (0:false - 1: true)
a = 0.03
b = -2
c = -50
d = 100
Iin = 0
alpha = 0.3
Vi = -41.0
Vmin = -37.0
Ka = 7
Ki = 8.75
tau = 6.5
celltype = 1 : A flag for indicating what kind of cell it is, used for changing the dynamics slightly (see list of cell types in initial comment).
alive = 1 : A flag for deciding whether or not the cell is alive -- if it's dead, acts normally except it doesn't fire spikes
cellid = -1 : A parameter for storing the cell ID, if required (useful for diagnostic information)
}
: Variables used for internal calculations
ASSIGNED {
v (mV)
i (nA)
u (mV) : Slow current/recovery variable
vt (mV): Dynamic threshold
delta
t0
derivtype
}
: Initial conditions
INITIAL {
u = u0
if (dyn_th == 1) {
vt = vt_thetainf(vr)
} else {
vt = vt0
}
derivtype=2
net_send(0,1) : Required for the WATCH statement to be active; v=vr initialization done there
}
: Define neuron dynamics
BREAKPOINT {
delta = t-t0 : Find time difference
if (celltype<5) {
if (t>0) {
u = u + delta*a*(b*(v-vr)-u) : Calculate recovery variable
if (dyn_th == 1) {
vt = vt + delta*(vt_thetainf(v) - vt)/tau
}
}
}
else {
: For FS neurons, include nonlinear U(v): U(v) = 0 when v<vb ; U(v) = 0.025(v-vb) when v>=vb (d=vb=-55)
if (celltype==5) {
if (v<d) {
u = u + delta*a*(0-u)
}
else {
u = u + delta*a*((0.025*(v-d)*(v-d)*(v-d))-u)
}
}
: For TC neurons, reset b
if (celltype==6) {
if (v>-65) {b=0}
else {b=15}
u = u + delta*a*(b*(v-vr)-u) : Calculate recovery variable
}
: For TRN neurons, reset b
if (celltype==7) {
if (v>-65) {b=2}
else {b=10}
u = u + delta*a*(b*(v-vr)-u) : Calculate recovery variable
}
}
t0=t : Reset last time so delta can be calculated in the next time step
i = -(k*(v-vr)*(v-vt) - u + Iin)/C/1000
}
FUNCTION vt_thetainf (v(mV)) (mV){
vt_thetainf = alpha * (v - Vi) + Vmin + Ka * log ( 1 + exp((v - Vi)/Ki) )
}
FUNCTION derivfunc () {
if (celltype==5 && derivtype==2) { : For FS neurons, include nonlinear U(v): U(v) = 0 when v<vb ; U(v) = 0.025(v-vb) when v>=vb (d=vb=-55)
derivfunc = a*(0-u)
} else if (celltype==5 && derivtype==1) { : For FS neurons, include nonlinear U(v): U(v) = 0 when v<vb ; U(v) = 0.025(v-vb) when v>=vb (d=vb=-55)
derivfunc = a*((0.025*(v-d)*(v-d)*(v-d))-u)
} else if (celltype==5) {
VERBATIM
hoc_execerror("izhi2007b.mod ERRA: derivtype not set",0);
ENDVERBATIM
} else {
derivfunc = a*(b*(v-vr)-u) : Calculate recovery variable
}
}
: Input received
NET_RECEIVE (w) {
: Check if spike occurred
if (flag == 1) { : Fake event from INITIAL block
if (celltype == 4) { : LTS cell
WATCH (v>(vpeak-0.1*u)) 2 : Check if threshold has been crossed, and if so, set flag=2
} else if (celltype == 6) { : TC cell
WATCH (v>(vpeak+0.1*u)) 2
} else { : default for all other types
WATCH (v>vpeak) 2
}
: additional WATCHfulness
if (celltype==6 || celltype==7) {
WATCH (v> -65) 3 : change b param
WATCH (v< -65) 4 : change b param
}
if (celltype==5) {
WATCH (v> d) 3 : going up
WATCH (v< d) 4 : coming down
}
v = v0 : initialization can be done here
: FLAG 2 Event created by WATCH statement -- threshold crossed for spiking
} else if (flag == 2) {
if (alive) {net_event(t)} : Send spike event if the cell is alive
: For LTS neurons
if (celltype == 4) {
v = c+0.04*u : Reset voltage
if ((u+d)<670) {u=u+d} : Reset recovery variable
else {u=670}
}
: For FS neurons (only update v)
else if (celltype == 5) {
v = c : Reset voltage
}
: For TC neurons (only update v)
else if (celltype == 6) {
v = c-0.1*u : Reset voltage
u = u+d : Reset recovery variable
} else {: For RS, IB and CH neurons, and RTN
v = c : Reset voltage
u = u+d : Reset recovery variable
}
: FLAG 3 Event created by WATCH statement -- v exceeding set point for param reset
} else if (flag == 3) {
: For TC neurons
if (celltype == 5) { derivtype = 1 : if (v>d) u'=a*((0.025*(v-d)*(v-d)*(v-d))-u)
} else if (celltype == 6) { b=0
} else if (celltype == 7) { b=2
}
: FLAG 4 Event created by WATCH statement -- v dropping below a setpoint for param reset
} else if (flag == 4) {
if (celltype == 5) { derivtype = 2 : if (v<d) u==a*(0-u)
} else if (celltype == 6) { b=15
} else if (celltype == 7) { b=10
}
}
}
