https://www.physiology.org/doi/full/10. ... b%3Dpubmed
The mechanism is governed by two equations: dv/dt for the membrane potential (v) and dtheta/dt for the firing threshold (theta).
Here is my code:
Code: Select all
TITLE esser_mech.mod
NEURON {
SUFFIX esser
POINT_PROCESS esser_mech
USEION na READ ena WRITE ina
USEION k READ ek WRITE ik
RANGE gNa_leak, gK_leak, theta_eq, tau_theta, tau_spike, tau_m, gspike, C, tspike, vrefrac, ifake
NONSPECIFIC_CURRENT ispike
}
UNITS {
(mV) = (millivolt)
(mA) = (milliamp)
(S) = (siemens)
}
PARAMETER {
gNa_leak = 0.14 (S/cm2)
gK_leak = 1.0 (S/cm2)
C = 0.85
theta_eq = -53 (mV) : resting threshold
tau_theta = 2 (ms) : threshold time constant
tau_spike = 1.75 (ms) : time const during a spike
tau_m = 15 (ms) : membrane time const
tspike = 2 (ms) : duration of spike
vap = 30 (mV) : v of spike
}
ASSIGNED {
v (mV)
ena (mV)
ek (mV)
ina (mA/cm2)
ik (mA/cm2)
ispike (mA/cm2)
gspike (S/cm2)
ifake (mA/cm2) : causes solver to execute the breakpoint block at the
: same time as the calculation of the other currents
}
STATE {
theta (mV)
}
BREAKPOINT {
SOLVE states METHOD cnexp
ina = gNa_leak*(v-ena)/tau_m : Sodium-like leak current
ik = gK_leak*(v-ek)/tau_m : Potassium-like leak current
ispike = gspike*(v-ek)/tau_spike
ifake = 0
}
INITIAL {
net_send(0,1)
theta = theta_eq
gspike = 0
ena = 30 (mV)
ek = -90 (mV)
v = -75.263 (mV)
}
DERIVATIVE states {
theta' = (-(theta - theta_eq) + C*(v - theta_eq))/tau_theta : threshold
}
NET_RECEIVE (w) {
if (flag == 1){
WATCH (v > theta) 2
}else if (flag == 2){
net_event(t)
gspike = 1
v = vap
theta = vap
net_send(tspike,3)
}else if (flag == 3){
gspike = 0
net_send(0,1)
}
}
My concerns are:
1. Is my implementation correct?
2. I know v is reserved for membrane potential, which is computed by numerical integration of the discretized cable equation. dv/dt equation in the paper (with a tau_m term) is not the same as the cable equation in NEURON. Would v be calculated correctly? Or do I need to define another variable for the membrane potential?
Thank you!